Coweeta Hydrologic Laboratory Watershed 32 using Static, SSURGO, and DSM soil inputs. Previously WS27.

The following R Markdown script prepares and runs the Regional Hydro-Ecologic Simulation System (RHESSYS) within the R environment. Modified from RHESSysIOinR example scripts created by Will Burke and Ryan Bart. Includes CLHS code snippets and SDA code from Dylan Beaudette.

Windows Subsystem via Linux must be installed before RHESSys will run in a Windows environment. GRASS8 must be installed to run the spatial preprocessing portion of the script.

Code has been modified to run RHESSys 7.4 and 5.18. 7.4 vegetation processing differs slightly from 5.18. Various template differences exist between the two versions of RHESSys and template files must be carefully reviewed if moving between versions. This script is currently set up to run RHESSys 7.4

#install.packages("devtools")
#devtools::install_github("ropensci/FedData")

library(aqp)
## This is aqp 1.41
library(daymetr)
library(clhs)
## Registered S3 methods overwritten by 'tibble':
##   method     from  
##   format.tbl pillar
##   print.tbl  pillar
library(EcoHydRology)
## Loading required package: operators
## 
## Attaching package: 'operators'
## The following objects are masked from 'package:base':
## 
##     options, strrep
## Loading required package: topmodel
## Loading required package: DEoptim
## Loading required package: parallel
## 
## DEoptim package
## Differential Evolution algorithm in R
## Authors: D. Ardia, K. Mullen, B. Peterson and J. Ulrich
## Loading required package: XML
library(elevatr)
library(ezknitr)
library(FedData)
## 
## Attaching package: 'FedData'
## The following object is masked from 'package:operators':
## 
##     %>%
library(foreach)
library(foreign)
library(ggplot2)
library(hydroGOF)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## 
## Attaching package: 'hydroGOF'
## The following object is masked from 'package:topmodel':
## 
##     NSeff
library(imputeTS)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## 
## Attaching package: 'imputeTS'
## The following object is masked from 'package:zoo':
## 
##     na.locf
## The following object is masked from 'package:operators':
## 
##     %>%
library(knitr)
library(latticeExtra)
## Loading required package: lattice
## 
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
## 
##     layer
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
library(Metrics)
## 
## Attaching package: 'Metrics'
## The following objects are masked from 'package:hydroGOF':
## 
##     mae, mse, rmse
library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following objects are masked from 'package:aqp':
## 
##     metadata, metadata<-
library(rasterVis)
library(readxl)
library(reshape2)
library(rgdal)
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.1, released 2021/12/27
## Path to GDAL shared files: C:/Users/Carlos/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/Carlos/Documents/R/win-library/4.0/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-6
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(rgrass)
## GRASS GIS interface loaded with GRASS version: (GRASS not running)
library(Rmisc)
## Loading required package: plyr
library(RHESSysPreprocessing)
library(RHESSysIOinR)
## 
## Attaching package: 'RHESSysIOinR'
## The following object is masked from 'package:RHESSysPreprocessing':
## 
##     read_world
library(sf)
## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE
## 
## Attaching package: 'sf'
## The following object is masked from 'package:operators':
## 
##     %>%
library(soilDB)
library(sp)
library(stringi)
library(tactile)
library(terra)
## terra 1.5.21
## 
## Attaching package: 'terra'
## The following object is masked from 'package:rgdal':
## 
##     project
## The following object is masked from 'package:knitr':
## 
##     spin
## The following object is masked from 'package:zoo':
## 
##     time<-
## The following object is masked from 'package:ggplot2':
## 
##     arrow
library(testthat)
## 
## Attaching package: 'testthat'
## The following object is masked from 'package:terra':
## 
##     describe
## The following object is masked from 'package:operators':
## 
##     %>%
library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v tibble  3.0.3     v dplyr   1.0.5
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.1.2     v forcats 0.5.1
## v purrr   0.3.4
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x forcats::%>%()           masks stringr::%>%(), dplyr::%>%(), purrr::%>%(), tidyr::%>%(), tibble::%>%(), testthat::%>%(), sf::%>%(), imputeTS::%>%(), FedData::%>%(), operators::%>%()
## x purrr::accumulate()      masks foreach::accumulate()
## x dplyr::arrange()         masks plyr::arrange()
## x terra::arrow()           masks ggplot2::arrow()
## x lubridate::as.difftime() masks base::as.difftime()
## x dplyr::combine()         masks aqp::combine()
## x purrr::compact()         masks plyr::compact()
## x dplyr::count()           masks plyr::count()
## x lubridate::date()        masks base::date()
## x readr::edition_get()     masks testthat::edition_get()
## x tidyr::extract()         masks terra::extract(), raster::extract()
## x dplyr::failwith()        masks plyr::failwith()
## x dplyr::filter()          masks stats::filter()
## x dplyr::id()              masks plyr::id()
## x terra::intersect()       masks raster::intersect(), lubridate::intersect(), base::intersect()
## x purrr::is_null()         masks testthat::is_null()
## x dplyr::lag()             masks stats::lag()
## x latticeExtra::layer()    masks ggplot2::layer()
## x readr::local_edition()   masks testthat::local_edition()
## x dplyr::matches()         masks tidyr::matches(), testthat::matches()
## x dplyr::mutate()          masks plyr::mutate()
## x dplyr::rename()          masks plyr::rename()
## x dplyr::select()          masks raster::select()
## x lubridate::setdiff()     masks base::setdiff()
## x dplyr::slice()           masks aqp::slice()
## x dplyr::src()             masks terra::src()
## x dplyr::summarise()       masks plyr::summarise()
## x dplyr::summarize()       masks plyr::summarize()
## x terra::union()           masks raster::union(), lubridate::union(), base::union()
## x purrr::when()            masks foreach::when()
library(viridisLite)


knitr::opts_knit$set(root.dir = system.file("extdata/", package = "RHESSysIOinR"))

## Prep RHESSys
expect_file_exists = function(path) {
  # 1. Capture object and label
  act <- quasi_label(rlang::enquo(path), arg = "path")
  # 2. Call expect()
  expect(
    file.exists(act$val),
    sprintf("%s file does not exist.", act$lab)
  )
  # 3. Invisibly return the value
  invisible(act$val)
}

expect_file_sizeKB_gt = function(path, size_KB) {
  # 1. Capture object and label
  act <- quasi_label(rlang::enquo(path), arg = "path")
  # 2. Call expect()
  expect(
    (file.size(act$val) / 1024) > size_KB,
    sprintf("%s file is not greater than %f KB.", act$lab, size_KB)
  )
  # 3. Invisibly return the value
  invisible(act$val)
}

setwd(system.file("extdata/", package = "RHESSysIOinR"))

### Downloads most recent version of RHESSYs
# gert::git_clone(url = "https://github.com/RHESSys/RHESSys", path = "./rh_dev", branch = "master")

### Make sure to change files in rh_dev to 5.18 if using older RHESSys files
# compile_rhessys(location = "rh_dev/")

rh_ver = dir(path = "rh_dev/rhessys/", pattern = "^rhessys\\d+",recursive = F)

test_that("compile_rhessys works + rhessys compiles via R system", {
  expect_gt(length(rh_ver), 0)
})
## Test passed
rh_path = file.path("rh_dev/rhessys/", rh_ver)

Load in Shapefiles for Watershed Setup

getwd()
## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"
setwd(system.file("extdata/", package = "RHESSysIOinR"))

watershedshapefile <-("C:/Users/Carlos/Desktop/RHESSys_R_Processing/Boundary Shapefiles/Coweeta_Hydrologic_Laboratory.shp")
streamshapefile<- ("C:/Users/Carlos/Desktop/ORISE/GIS/CW/Streamflow/StreamsCaldwellClipped.shp")
weirshapefile<- ("C:/Users/Carlos/Desktop/ORISE/GIS/CW/Weirs/Subbasinoutlets.shp")
roadshapefile <- ( "C:/Users/Carlos/Desktop/ORISE/GIS/CW/Roads/UpdatedRoads4Reprojected32617.shp")

#streamflow dataset
#soil moisture dataset
##Quick way to create base file necessary for climate processing not functional yet

Download Daymet Data for Watershed of Interest using daymetr package Currently only creates precip, tmin and tmax files.

Extend DAYMET datasets to allow for longer spinup period Read Daymet Datasets and Replicate 5x times to create a longer period of record for spin up Create base File manually for now

## this may be the source of lag in outputs/measured data

getwd()
## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"
precip <- read.csv('clim/cwtws32daymet.rain', header = FALSE)
tmin <- read.csv('clim/cwtws32daymet.tmin', header = FALSE)
tmax <- read.csv('clim/cwtws32daymet.tmax', header = FALSE)
#tavg <- read.csv('clim/cwtws32daymet.tavg', header = FALSE)
#rh <- read.csv('clim/cwtws32daymet.relative_humidity', header = FALSE)
#wind <- read.csv('clim/cwtws32daymet.wind', header = FALSE)

head(precip)
##           V1
## 1 1980 1 1 1
## 2    0.00427
## 3          0
## 4    0.00866
## 5    0.01468
## 6          0
nrow(precip)-1
## [1] 15341
substr(precip[1,1],1,nchar(precip[1,1])-2)
## [1] "1980 1 1"
dt <- as.Date(substr(precip[1,1],1,nchar(precip[1,1])-2), "%Y %m %d")

dt.new<- dt - as.difftime((nrow(precip)-1)*5, units="days")

#remove first value and repeats dataset 6 times

longcwtraindata<- do.call("rbind",replicate(6,as.list((precip[-1,]))))
longcwttmindata<- do.call("rbind",replicate(6,as.list((tmin[-1,]))))
longcwttmaxdata<- do.call("rbind",replicate(6,as.list((tmax[-1,]))))
#longcwttavgdata<- do.call("rbind",replicate(6,as.list((tavg[-1,]))))
#longcwtrhdata<- do.call("rbind",replicate(6,as.list((rh[-1,]))))
#longcwtwinddata<- do.call("rbind",replicate(6,as.list((wind[-1,]))))

longcwstartdate<- paste(year(ymd(dt.new)),month(ymd(dt.new)),day(ymd(dt.new)),"1")

longcwtrain<- rbind(longcwstartdate,longcwtraindata)
longcwttmin<- rbind(longcwstartdate,longcwttmindata)
longcwttmax<- rbind(longcwstartdate,longcwttmaxdata)
#longcwttavg<- rbind(longcwstartdate,longcwttavgdata)
#longcwtrh<- rbind(longcwstartdate,longcwtrhdata)
#longcwtwind<- rbind(longcwstartdate,longcwtwinddata)

longcwtrain<- data.frame(longcwtrain)
longcwttmin<- data.frame(longcwttmin)
longcwtmax<- data.frame(longcwttmax)
#longcwttavg<- data.frame(longcwttavg)
#longcwtrh<- data.frame(longcwtrh)
#longcwtwind<- data.frame(longcwtwind)

write.table(longcwtrain, file = "clim/cwtws32extended.rain", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(longcwttmin, file = "clim/cwtws32extended.tmin", quote = FALSE, row.names = FALSE, col.names = FALSE)
write.table(longcwttmax, file = "clim/cwtws32extended.tmax", quote = FALSE, row.names = FALSE, col.names = FALSE)
#write.table(longcwttavg, file = "clim/cwtws32extended.tavg", quote = FALSE, row.names = FALSE, col.names = FALSE)
#write.table(longcwtrh, file = "clim/cwtws32extended.relative_humidity", quote = FALSE, row.names = FALSE, col.names = FALSE)
#write.table(longcwtwind, file = "clim/cwtws32extended.wind", quote = FALSE, row.names = FALSE, col.names = FALSE)

read_clim("clim/cwtws32extended.rain", dates_out=TRUE)
## [1] "1769 12 27 01" "2021 12 31 24"
read_clim("clim/cwtws32daymet.rain", dates_out=TRUE)
## [1] "1980 01 01 01" "2021 12 31 24"
#daymet_data$data

# precip data for watershed not readily available

Download and Prepare DEM Code to download DEM originally prepared by Dylan Beaudette

# use watershed polygons for BBOX to request elevation data
x <- read_sf(watershedshapefile)
# buffer to ensure that there are no truncated watersheds
x.buff <- st_as_sf(st_buffer(st_as_sfc(st_bbox(x)), dist = 1000))
# convert to GCS WGS84 -> no transformation / warp will be applied to DEM
x.gcs <- st_transform(x, 4326)
x.buff.gcs <- st_transform(x.buff, 4326)
# requires sf collection (geometry + attributes)
# get DEM in GCS WGS84 and use z = 14 for best available data
e <- get_elev_raster(locations = x.buff.gcs, z = 14, clip = 'bbox')
## Mosaicing & Projecting
## Clipping DEM to bbox
## Note: Elevation units are in meters.
res(e)
## [1] 3.921794e-05 3.921794e-05
# check: OK
{plot(e, main = "DEM, Watershed Outline and Buffer" )
  plot(st_geometry(x.buff.gcs), add = TRUE, col = NA, lwd = 2)
  plot(st_geometry(x.gcs), add = TRUE, col = NA)}

{plot(e, main = "DEM and selected Watershed")
  plot(st_geometry(x.gcs[which(x$WS=='32'),]), add = TRUE, col = NA)}

cropped_dem<- crop(e,st_bbox(x.gcs))

plot(cropped_dem, main = "DEM cropped to Buffer")
plot(st_geometry(x.gcs[which(x$WS=='32'),]), add = TRUE, col = NA)

newproj<- "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs"

reproje<- projectRaster(from = cropped_dem, crs = newproj, method = 'ngb')

plot(reproje, main = "Reprojected cropped DEM")

dummydata<- reproje

res(dummydata) <- 10

resampleddem <- resample(reproje, dummydata, method = 'bilinear')

plot(resampleddem, main = "DEM reprojected to 10m Resolution")

Start Spatial Processing in GRASS8 Code to prepare spatial input files for RHESSYs takes cues from Will Burkes “spatial_input_gen.R” script but is prepared for GRASS8 instead of GRASS7

#gisbase_grass <- ifelse(.Platform$OS.type == "windows", link2GI::paramGRASSw()$gisbase_GRASS[1], link2GI::paramGRASSx()$gisbase_GRASS[1])

## set easting and northing as coordinates of lowest point within boundary basin, set threshold for basin delineation
threshold = 400
easting = NULL
northing = NULL
 

initGRASS(gisBase = "C:/Program Files/GRASS GIS 8.2",
          gisDbase = "C:/Users/Carlos/Documents/grassdata",
          location = "Coweeta", 
          mapset = "PERMANENT", 
          override = TRUE)
## Warning in system(syscmd, intern = intern, ignore.stderr = ignore.stderr, :
## running command 'g.proj.exe -w' had status 884
## gisdbase    C:/Users/Carlos/Documents/grassdata 
## location    Coweeta 
## mapset      PERMANENT 
## rows        611 
## columns     537 
## north       3884060 
## south       3877950 
## west        273846.7 
## east        279216.7 
## nsres       10 
## ewres       10 
## projection  +proj=utm +no_defs +zone=17 +a=6378137 +rf=298.257223563
## +towgs84=0.000,0.000,0.000 +type=crs +to_meter=1
execGRASS("g.proj", proj4="+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs", flag = "c" )

execGRASS("g.gisenv", flag="n")

execGRASS("g.region", flag=c("p"))

##Convert DEM to Spatrast and then import DEM to GRASS 

spatrast_dem10 <- rast(resampleddem)

write_RAST(spatrast_dem10, "grass_dem10", overwrite = TRUE)
## SpatRaster read into GRASS using r.in.gdal from memory
## The following code is heavily influenced by Will Burkes "spatial_input_gen.R" script

execGRASS("g.region", raster = "grass_dem10", flag=c("p"))

## Import shapefiles (Streams, Weirs, Watersheds, Roads)
## only roads and Watersheds necessary to Run

streams<- vect("C:/Users/Carlos/Desktop/ORISE/GIS/CW/Streamflow/StreamsCaldwellClipped.shp")
streams<- project(streams,"+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs")
write_VECT(streams, "streams")

weirs<- vect("C:/Users/Carlos/Desktop/ORISE/GIS/CW/Weirs/Subbasinoutlets.shp")
weirs<- project(weirs,"+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs")
write_VECT(weirs, "weirs")

watersheds<- vect("C:/Users/Carlos/Desktop/RHESSys_R_Processing/Boundary Shapefiles/Coweeta_Hydrologic_Laboratory.shp")
wgs84watersheds<- project(watersheds,"+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs")
write_VECT(wgs84watersheds, "watersheds")

roads <- vect( "C:/Users/Carlos/Desktop/ORISE/GIS/CW/Roads/UpdatedRoads4Reprojected32617.shp")
## Warning: [vect] Z coordinates ignored
roads<- project(roads,"+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs")
write_VECT(roads,"roads")

## Check all vectors that have been exported to GRASS

execGRASS("g.list", parameters = list(type = "vector"))

## Generate a depressionless DEM and d-8 flow direction maps from DEM

execGRASS("r.fill.dir", input="grass_dem10", output="dem10f", direction="dir10", flags = c("overwrite") )

## Generate flow accumulation, drainage, and horizon maps

execGRASS("r.watershed", elevation = "dem10f", accumulation = "acc10", drainage = "drain10")

## Generate slope, aspect, wetness index, horizon, and road maps

execGRASS("r.slope.aspect", elevation = "grass_dem10", slope = "slope10", aspect = "aspect10")

execGRASS("r.topidx", input = "grass_dem10", output = "topidx10")

execGRASS("r.mapcalc", expression = "topidx10_100 = topidx10 * 100", flags = "overwrite")

execGRASS("r.horizon", flags = "d", elevation = "grass_dem10", direction = 0, output = "east")
execGRASS("r.horizon", flags = "d", elevation = "grass_dem10", direction = 180, output = "west")

execGRASS("r.mapcalc", expression = "e_horizon_1000 = sin(east_000) * 1000", flags = "overwrite")
execGRASS("r.mapcalc", expression = "w_horizon_1000 = sin(west_180) * 1000", flags = "overwrite")

## Generate a hillslope and stream map

subbasin_name = paste0("subbasin_th",as.character(threshold))
stream_name = paste0("stream_th",as.character(threshold))
hillslope_name = paste0("hillslope_th",as.character(threshold))

execGRASS("r.watershed", elevation = "dem10f", threshold = threshold, basin = subbasin_name, stream = stream_name, half_basin = hillslope_name)

## Delineate watershed boundary

weirs[12,]
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 1, 1  (geometries, attributes)
##  extent      : 276120.4, 276120.4, 3881339, 3881339  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=utm +zone=17 +datum=WGS84 +units=m +no_defs 
##  names       :    id
##  type        : <int>
##  values      :    12
easting = 276122.49
northing = 3881345.46

execGRASS("g.region", raster = "grass_dem10", flag=c("p"))

# r.water.outlet input=drain10 output=basin_ws321 coordinates=276122.49,3881345.46 --overwrite

execGRASS("r.water.outlet", input ="drain10", output="basin_ws32",coordinates = c(easting, northing), flags = "overwrite")

## Mask Data
# execGRASS("r.mask", input = "basin_ws32", maskcats = 1)

execGRASS("r.info",map = "basin_ws32")

execGRASS("r.info",map = "grass_dem10")

## Generate Patch Maps

execGRASS("r.mapcalc", expression = "patch = row() * 537 + col()", flags = "overwrite")


## Create Road Maps

execGRASS("v.to.rast", input="roads" ,output="roads", use="cat", label_column="ROADNAME", flags = "overwrite")

execGRASS("r.mapcalc", expression = "roads = roads > 0", flags = "overwrite")

execGRASS("r.null", map = "roads", null= 0)

## Check all rasters that have been exported to or created in GRASS

execGRASS("g.list", parameters = list(type = "raster"))

## Export Prepared Maps from GRASS to R

filleddem <- read_RAST("dem10f")

##Export maps back into R or export directly from GRASS to disk

grass_acc10 <- read_RAST("acc10")
grass_aspect10 <- read_RAST("aspect10")
grass_basin_ws32 <- read_RAST("basin_ws32")
grass_dem10 <- read_RAST("grass_dem10")
grass_dem10f <- read_RAST("dem10f")
grass_dir10 <- read_RAST("dir10")
grass_drain10 <- read_RAST("drain10")
grass_e_horizon_1000 <- read_RAST("e_horizon_1000")
grass_hillslope_th1000 <- read_RAST("hillslope_th1000")
#grass_patch <- read_RAST("patch")
grass_slope10 <- read_RAST("slope10")
grass_stream_th1000 <- read_RAST("stream_th1000")
grass_topidx10 <- read_RAST("topidx10")
grass_topidx10_100 <- read_RAST("topidx10_100")
grass_w_horizon_1000 <- read_RAST("w_horizon_1000")
grass_east_000 <- read_RAST("east_000")
grass_west_180 <- read_RAST("west_180")
grass_subbasin_th1000 <- read_RAST("subbasin_th1000")


execGRASS("r.out.gdal", input = "acc10", output = "grassexport/cwt/acc10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "aspect10", output = "grassexport/cwt/aspect10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "basin_ws32", output = "grassexport/cwt/basin_ws32.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "grass_dem10", output = "grassexport/cwt/grass_dem10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "dem10f", output = "grassexport/cwt/dem10f.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "dir10", output = "grassexport/cwt/dir10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "drain10", output = "grassexport/cwt/drain10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "e_horizon_1000", output = "grassexport/cwt/e_horizon_1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "hillslope_th1000", output = "grassexport/cwt/hillslope_th1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "patch", output = "grassexport/cwt/patch.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "slope10", output = "grassexport/cwt/slope10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "stream_th1000", output = "grassexport/cwt/stream_th1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "topidx10", output = "grassexport/cwt/topidx10.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "topidx10_100", output = "grassexport/cwt/topidx10_100.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "w_horizon_1000", output = "grassexport/cwt/w_horizon_1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "east_000", output = "grassexport/cwt/east_000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "west_180", output = "grassexport/cwt/west_180.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "subbasin_th1000", output = "grassexport/cwt/subbasin_th1000.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "roads", output = "grassexport/cwt/roads.tif", format = "GTiff", flags = "overwrite")

execGRASS("r.out.gdal", input = "stream_th1500", output = "grassexport/cwt/stream_th1500.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "stream_th500", output = "grassexport/cwt/stream_th500.tif", format = "GTiff", flags = "overwrite")
execGRASS("r.out.gdal", input = "stream_th400", output = "grassexport/cwt/stream_th400.tif", format = "GTiff", flags = "overwrite")

Crop Maps that were prepared in GRASS. Maps are cropped R because MASK command is not functioning properly for GRASS via R. Original input maps are been cropped before processing worldfile and flowtable in R, without cropping the RHESSYs worldfile and flowtable preparation in the R environment will not function properly,

ws32<- raster("grassexport/cwt/basin_ws32.tif")
acc10<- raster("grassexport/cwt/acc10.tif")
patches<- raster("grassexport/cwt/patch.tif")
aspect10<- raster("grassexport/cwt/aspect10.tif")
dem10<- raster("grassexport/cwt/grass_dem10.tif")
dem10f<- raster("grassexport/cwt/dem10f.tif")
dir10<- raster("grassexport/cwt/dir10.tif")
drain10<- raster("grassexport/cwt/drain10.tif")
hillslope<- raster("grassexport/cwt/hillslope_th1000.tif")
roads<- brick("grassexport/cwt/roads.tif")
slope10<- raster("grassexport/cwt/slope10.tif")
streams<- raster("grassexport/cwt/stream_th1000.tif")
topidx10_100<- raster("grassexport/cwt/topidx10_100.tif")
e_horizon_1000 <- raster("grassexport/cwt/e_horizon_1000.tif")
w_horizon_1000 <- raster("grassexport/cwt/w_horizon_1000.tif")
subbasins<- raster("grassexport/cwt/subbasin_th1000.tif")

ws32polygon<- rasterToPolygons(ws32, fun = function(x){x==1}, dissolve = TRUE, digits = 12)
## Loading required namespace: rgeos
maskedws32 <- crop(mask(ws32,ws32polygon),ws32polygon)
maskedacc10 <- crop(mask(acc10,ws32polygon),ws32polygon)
maskedpatches <- crop(mask(patches,ws32polygon),ws32polygon)
maskedaspect10 <- crop(mask(aspect10,ws32polygon),ws32polygon)
maskeddem10 <- crop(mask(dem10,ws32polygon),ws32polygon)
maskeddem10f <- crop(mask(dem10f,ws32polygon),ws32polygon)
maskeddir10 <- crop(mask(dir10,ws32polygon),ws32polygon)
maskeddrain10 <- crop(mask(drain10,ws32polygon),ws32polygon)
maskedhillslope <- crop(mask(hillslope,ws32polygon),ws32polygon)
maskedroads <- crop(mask(roads,ws32polygon),ws32polygon)
maskedslope10 <- crop(mask(slope10,ws32polygon),ws32polygon)
maskedstreams <- crop(mask(streams,ws32polygon),ws32polygon)
maskedtopidx10_100 <- crop(mask(topidx10_100,ws32polygon),ws32polygon)
maskede_horizon_1000 <- crop(mask(e_horizon_1000,ws32polygon),ws32polygon)
maskedw_horizon_1000 <- crop(mask(w_horizon_1000,ws32polygon),ws32polygon)
maskedsubbasins<- crop(mask(subbasins,ws32polygon),ws32polygon)

## Save as tif

writeRaster(maskedws32,"grassexport/cwt_ws32/basin_ws32",format = "GTiff", overwrite = TRUE)
writeRaster(maskedacc10,"grassexport/cwt_ws32/acc10",format = "GTiff", overwrite = TRUE)
writeRaster(maskedpatches,"grassexport/cwt_ws32/patch",format = "GTiff", overwrite = TRUE)
writeRaster(maskedaspect10,"grassexport/cwt_ws32/aspect10",format = "GTiff", overwrite = TRUE)
writeRaster(maskeddem10,"grassexport/cwt_ws32/dem10",format = "GTiff", overwrite = TRUE)
writeRaster(maskeddem10f,"grassexport/cwt_ws32/dem10f",format = "GTiff", overwrite = TRUE)
writeRaster(maskeddir10,"grassexport/cwt_ws32/dir10",format = "GTiff", overwrite = TRUE)
writeRaster(maskeddrain10,"grassexport/cwt_ws32/drain10",format = "GTiff", overwrite = TRUE)
writeRaster(maskedhillslope,"grassexport/cwt_ws32/hillslope",format = "GTiff", overwrite = TRUE)
writeRaster(maskedroads,"grassexport/cwt_ws32/roads",format = "GTiff", overwrite = TRUE)
writeRaster(maskedslope10,"grassexport/cwt_ws32/slope10",format = "GTiff", overwrite = TRUE)
writeRaster(maskedstreams,"grassexport/cwt_ws32/streams",format = "GTiff", overwrite = TRUE)
writeRaster(maskedtopidx10_100,"grassexport/cwt_ws32/topidx10_100",format = "GTiff", overwrite = TRUE)
writeRaster(maskede_horizon_1000,"grassexport/cwt_ws32/e_horizon_1000",format = "GTiff", overwrite = TRUE)
writeRaster(maskedw_horizon_1000,"grassexport/cwt_ws32/w_horizon_1000",format = "GTiff", overwrite = TRUE)
writeRaster(maskedsubbasins,"grassexport/cwt_ws32/subbasins",format = "GTiff", overwrite = TRUE)


plot(maskedws32, main ="ws32")

plot(maskedpatches, main = "patches")

plot(maskedacc10, main = "acc10")

plot(maskedaspect10, main = "aspect10")

plot(maskedsubbasins, main = "subbasins")

plot(maskeddem10, main = "dem10")

plot(maskeddem10f, main = "dem10f")

plot(maskeddir10, main = "dir10")

plot(maskeddrain10, main = "drain10")

plot(maskedhillslope, main = "hillslope")

plot(maskedroads, main = "roads")

plot(maskedslope10, main = "slope10")

plot(maskedstreams, main = "streams")

plot(maskedtopidx10_100, main = "topidx10_100")

plot(maskede_horizon_1000, main = "e_horizon_1000")

plot(maskedw_horizon_1000, main = "w_horizon_1000")

SSURGO-via-SDA Code originally prepared by Dylan Beaudette Download SSURGO data and make sure CRS/Extent/Crop are Correct for input into RHESSys

# load BBOX
x <- read_sf(watershedshapefile)
# get gSSURGO grid here
mu <- mukey.wcs(aoi = x, db = 'gssurgo')
# unique map unit keys
ll <- levels(mu)[[1]]
# note SSA boundary
levelplot(mu, att = 'ID', margin = FALSE, colorkey = FALSE, col.regions = viridis, main = "SSURGO Map Units")

## Note: some of these map units are dominated by non-soil components
# Dominant Component (Numeric) -> NODATA
# Weighted Average -> use any available data, but wt. mean are diluted (see source data)

# get thematic data from SDA
# dominant component
# depth-weighted average
# sand, silt, clay, AWC (RV)

p <-  get_SDA_property(property = c("sandtotal_r","silttotal_r","claytotal_r", "awc_r"),
                       method = "Weighted Average", 
                       mukeys = ll$ID,
                       top_depth = 0,
                       bottom_depth = 200)

head(p)
##   areasymbol musym
## 1      NC113   BuD
## 2      NC113   BuF
## 3      NC113   CaF
## 4      NC113   CdD
## 5      NC113   CdE
## 6      NC113   CdF
##                                                                           muname
## 1 Burton-Craggey-Rock outcrop complex, windswept, 15 to 30 percent slopes, stony
## 2 Burton-Craggey-Rock outcrop complex, windswept, 30 to 95 percent slopes, stony
## 3                     Cashiers gravelly fine sandy loam, 50 to 95 percent slopes
## 4                     Chandler gravelly fine sandy loam, 15 to 30 percent slopes
## 5                     Chandler gravelly fine sandy loam, 30 to 50 percent slopes
## 6                     Chandler gravelly fine sandy loam, 50 to 95 percent slopes
##    mukey sandtotal_r silttotal_r claytotal_r awc_r
## 1 545800       66.36       19.50       14.14  0.05
## 2 545801       66.41       19.28       14.32  0.05
## 3 545803       66.82       21.68       11.50  0.14
## 4 545805       46.74       41.76       11.50  0.13
## 5 545806       46.74       41.76       11.50  0.13
## 6 545807       46.74       41.76       11.50  0.13
# re-create raster attribute table with aggregate soil properties
rat <- merge(ll, p, by.x = 'ID', by.y = 'mukey', sort = FALSE, all.x = TRUE)
# re-pack RAT
levels(mu) <- rat
# convert raster + RAT --> stack of values
s <- deratify(mu, att = c("sandtotal_r", "silttotal_r", "claytotal_r", "awc_r"))
# graphical check
levelplot(
  s[[1:3]], 
  main = 'Sand, Silt, Clay (RV) 0-200cm\nWeighted Average',
  margin = FALSE, 
  scales = list(draw = FALSE), 
  col.regions = viridis
)

levelplot(
  s[[4]], 
  main = 'AWC (RV) 0-200cm\nWeighted Average',
  margin = FALSE, 
  scales = list(draw = FALSE), 
  col.regions = viridis
)

# convert to a representative soil texture class
txt.lut <- read.csv('http://soilmap2-1.lawr.ucdavis.edu/800m_grids/RAT/texture_2550.csv')
# make a copy
texture_060 <- s[[1]]
# note: soil textures that aren't present are dropped from factor levels
texture_060[] <- ssc_to_texcl(sand = s$sandtotal_r[], clay = s$claytotal_r[])
# extract RAT
rat <- levels(texture_060)[[1]]
# add colors
rat <- merge(rat, txt.lut, by.x = 'VALUE', by.y = 'class', all.x = TRUE, sort = FALSE)
# fix column order
rat <- rat[, c('ID', 'VALUE', 'hex', 'names')]
# re-pack
levels(texture_060) <- rat
# check: ok
cols <- levels(texture_060)[[1]]$hex
levelplot(texture_060, att = 'names', margin = FALSE, col.regions = cols, main = "Soil Textures") 

writeRaster(texture_060, file = 'C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt/ws32/soiltexture060.tif', overwrite = TRUE, options=c("COMPRESS=LZW"))

texture60raster<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt/ws32/soiltexture060.tif")

### this has been modified to use new SSURGO Inputs with deep-shallow depth classes

texture60raster<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt/ws32/ssurgo-soiltype-class.tif")

## reproject to raster mask instead

ws32<- raster("grassexport/cwt_ws32/ws32.tif")

newproj<- "+proj=utm +zone=17 +datum=WGS84 +units=m +no_defs"
reproj60 <- projectRaster(from = texture60raster, crs = proj4string(ws32), method = 'ngb')

plot(texture60raster, zlim = c(1,5), main = "Texture Raster with modified Depth Classes")

{plot(reproj60, zlim = c(1,5), main = "Reprojected Texture Raster with modified Depth Classes and Watersheds")
plot(st_geometry(x), add = TRUE, col = NA)}

## prepare 10m resolution raster and resample prepared soil data to 10m
resample10<- raster(resolution=c(10,10), crs = proj4string(reproj60), ext=extent(reproj60))

reproj6010 <- resample(reproj60,resample10, method = 'ngb')

{plot(reproj6010, zlim = c(1,5), main = "Reprojected and Resampled Texture Raster with modified Depth Classes and Watersheds")
plot(st_geometry(x), add = TRUE, col = NA)}

## Crop prepared raster to extent of watershed of interest

## croppedtexture60 <- crop(reproj6010,extent(x[which(x$WS=='32'),]))

#for now cropping is completed based on delineated watershed extent from raster instead of watershed extent from watershed polygons
croppedtexture60 <- crop(reproj6010,extent(ws32))

# croppedtexture60 <- projectRaster(from = croppedtexture60, crs = newproj, method = 'ngb')

## Mask prepared raster using watershed polygon
#masked60 <- mask(croppedtexture60,x[which(x$WS=='32'),])

## for now masking is completed based on delineated watershed extent from raster instead of watershed extent from watershed polygons
ws32polygon<- rasterToPolygons(ws32, fun = function(x){x==1}, dissolve = TRUE, digits = 12)


extent(croppedtexture60)<- extent(ws32)



{plot(reproj6010, zlim = c(1,5), main = "Reprojected and Resampled Texture Raster with modified Depth Classes and Selected Watershed")
plot(ws32polygon, add = TRUE, col = NA)}

masked60 <- mask(croppedtexture60,ws32polygon)

plot(masked60, zlim = c(1,5), main = "Reprojected and Resampled Texture Raster with modified Depth Classes and Masked Selected Watershed")

plot(x[which(x$WS=='32'),])

plot(st_geometry(x[which(x$WS=='32'),]))

{plot(reproj60, zlim = c(1,5))
plot(st_geometry(x[which(x$WS=='32'),]), add = TRUE, col = NA)}

{plot(as.factor(croppedtexture60), zlim = c(1,5))
plot(st_geometry(x[which(x$WS=='32'),]), add = TRUE, col = NA)}

{plot(as.factor(masked60), zlim = c(1,5))
plot(st_geometry(x[which(x$WS=='32'),]), add = TRUE, col = NA)}

{plot(masked60, zlim = c(1,5), main = "Polygon Delineation")
plot(st_geometry(x[which(x$WS=='32'),]), add = TRUE, col = NA)}

{plot(masked60, zlim = c(1,5), main = "Watershed Delineation")
plot(ws32polygon, add = TRUE, col = NA)}

Code Prepared By Dylan Beaudette and modified

setwd(system.file("extdata/", package = "RHESSysIOinR"))

ws32<- raster("grassexport/cwt_ws32/ws32.tif")

# RSS: cell values are map unit keys
r <- raster('Raster soil survey/Coweeta_Final_raster.tif')

## for now, convert to UTM z17
r <- projectRaster(r, ws32, method = 'ngb')

# convert to raster + RAT
r <- ratify(r)

# RAT: raster attribute table
rat <- read.dbf('Raster soil survey/Coweeta_Final_raster.tif.vat.dbf', as.is = TRUE)
names(rat) <- tolower(names(rat))

# musym is the national mu symbol
rss.mu <- st_read('Raster soil survey/rss_nc/rss_nc.gdb', layer = 'mapunit')
## Reading layer `mapunit' from data source 
##   `C:\Users\Carlos\Documents\R\win-library\4.0\RHESSysIOinR\extdata\Raster soil survey\rss_nc\rss_nc.gdb' 
##   using driver `OpenFileGDB'
rss.co <- st_read('Raster soil survey/rss_nc/rss_nc.gdb', layer = 'component')
## Reading layer `component' from data source 
##   `C:\Users\Carlos\Documents\R\win-library\4.0\RHESSysIOinR\extdata\Raster soil survey\rss_nc\rss_nc.gdb' 
##   using driver `OpenFileGDB'
rss.hz <- st_read('Raster soil survey/rss_nc/rss_nc.gdb', layer = 'chorizon')
## Reading layer `chorizon' from data source 
##   `C:\Users\Carlos\Documents\R\win-library\4.0\RHESSysIOinR\extdata\Raster soil survey\rss_nc\rss_nc.gdb' 
##   using driver `OpenFileGDB'
## check for missing symbols
# none missing from RAT
setdiff(rat$MUSYM, rss.mu$musym)
## NULL
# map unit table contains some extra
setdiff(rss.mu$musym, rat$MUSYM)
##  [1] "2y9mf" "2y9mh" "2y9mj" "2y9mk" "2y9ml" "2y9mm" "2y9mn" "2y9mp" "2y9mq"
## [10] "2y9mr" "2y9ms" "2y9mt" "2y9mv" "2y9mw" "2y9mx" "2y9mg" "2y9pb" "2y9pc"
## [19] "2y9pd" "2y9p9" "2y9nz" "2y9nl" "2y9nm" "2y9nw" "2y9nx" "2zyxq" "2zyxw"
## [28] "2zyxy" "2zyy5" "2zyxz" "2zyy0" "2zyxx" "2zyxt" "2zyy6" "2y9np" "2y9nt"
## [37] "2y9nv" "2zyxr" "2zyy1" "2zyy2" "2zyy3" "2y9nh" "2y9nj" "2y9nk" "2zyxv"
## [46] "2zyyb" "2zyyc" "2zyyd" "2zyyf" "2y9nn" "2y9p1" "2y9p2" "2y9nq" "2y9nr"
## [55] "2zyy4"
## subset child tables
rss.mu <- rss.mu[rss.mu$musym %in% rat$musym, ]
rss.co <- rss.co[rss.co$mukey %in% rss.mu$mukey, ]
rss.hz <- rss.hz[rss.hz$cokey %in% rss.co$cokey, ]

.dominantCondition <- function(i, v) {
  
  i <- i[i$compkind != 'Miscellaneous area', ]
  if(nrow(i) < 1) {
    return(NULL)
  }
  
  fm <- as.formula(sprintf("comppct_r ~ %s", v))
  a <- aggregate(fm, data = i, FUN = sum, na.rm = TRUE)
  
  idx <- order(a[['comppct_r']], decreasing = TRUE)[1]
  
  res <- data.frame(
    mukey = i$mukey[1],
    v = a[[v]][idx]
  )
  names(res) <- c('mukey', v)
  
  return(res)
}

.dominantValue <- function(i, v) {
  
  i <- i[i$compkind != 'Miscellaneous area', ]
  if(nrow(i) < 1) {
    return(NULL)
  }
  
  
  idx <- order(i[['comppct_r']], decreasing = TRUE)[1]
  
  res <- data.frame(
    mukey = i$mukey[1],
    v = i[[v]][idx]
  )
  
  names(res) <- c('mukey', v)
  
  return(res)
  
}

getDominantCondition <- function(x, v) {
  
  res <- lapply(
    split(rss.co, rss.co$mukey), 
    .dominantCondition, v = v
  )
  
  res <- do.call('rbind', res)
  
  return(res)
}

getDominantValue <- function(x, v) {
  
  res <- lapply(
    split(rss.co.aws050, rss.co.aws050$mukey), 
    .dominantValue, v = v
  )
  
  res <- do.call('rbind', res)
  
  return(res)
}

co.taxpartsize <- getDominantCondition(rss.co, v = 'taxpartsize')


co.spc <- rss.hz
depths(co.spc) <- cokey ~ hzdept_r + hzdepb_r
hzdesgnname(co.spc) <- 'hzname'

# at dZ = 1cm, use awc_r directly
co.spc$aws <- co.spc$awc_r * 1

a <- slab(co.spc, cokey ~ aws, slab.structure = c(0, 50), slab.fun = sum)
## Note: aqp::slice() will be deprecated in aqp version 2.0
## --> Please consider using the more efficient aqp::dice()
head(a)
##   variable           cokey value top bottom contributing_fraction
## 1      aws 3244721:2807027  2.40   0     50                     1
## 2      aws 3244721:2807439  4.02   0     50                     1
## 3      aws 3244721:2807442  2.50   0     50                     1
## 4      aws 3244721:2807445  3.73   0     50                     1
## 5      aws 3244722:2760595  2.40   0     50                     1
## 6      aws 3244722:2760596  4.02   0     50                     1
rss.co.aws050 <- merge(rss.co, a, by = 'cokey', sort = FALSE)[, c('mukey', 'cokey', 'compkind', 'comppct_r', 'value')]

# no NA allowed in wt. mean / etc.
rss.co.aws050 <- rss.co.aws050[!is.na(rss.co.aws050$value), ]

co.aws050 <- getDominantValue(rss.co.aws050, v = 'value')

names(co.aws050) <- c('mukey', 'aws050')


mu.subset <- rss.mu[, c('mukey', 'musym', 'mukind')]

agg <- merge(mu.subset, co.taxpartsize, by = 'mukey', sort = FALSE)

agg <- merge(agg, co.aws050, by = 'mukey', sort = FALSE)

rat <- merge(rat, agg, by = 'musym', sort = FALSE)

# attempt to split out component / series names
rat$co.names <- stri_split_fixed(str = rat$muname,  pattern = ',', n = 2, simplify = TRUE)[, 1]
# this only works because all component names are single-word names
rat$co.names <- stri_split_fixed(str = rat$co.names,  pattern = ' ', n = 2, simplify = TRUE)[, 1]
# RAT management
ll.original <- levels(r)[[1]]
# merge and pack updated RAT
ll <- merge(ll.original, rat, by.x = 'ID', by.y = 'value', all.x = TRUE, sort = FALSE)
levels(r) <- ll

# convert select attributes to new raster objects via RAT + codes

r.taxpartsize <- deratify(r, att = 'taxpartsize')


plot(r.taxpartsize)

ws32polygon<- rasterToPolygons(ws32, fun = function(x){x==1}, dissolve = TRUE, digits = 12)

maskeddsm <- mask(r.taxpartsize,ws32polygon)

plot(maskeddsm, main = "Masked DSM")

Reclassify Soil Maps

Default Texture Values in RHESSys 1 - Clay 2 - Silty Clay 3 - Silty Clay Loam 4 - Sandy Clay 5 - Sandy Clay Loam 6 - Clay Loam 7 - Silt 8 - Silty Loam 9 - Loam 10 - Sand 11 - Loamy Sand 12 - Sandy Loam 13 - Rock 14 - Water

Shallow Files Created and added to defs folder 19 - Shallow Sandy Clay Loam 23 - Shallow Loam 26 - Shallow Sandy Loam

Default Texture Values used by SDA 1 - loamy sand 2 - sandy loam 3 - loam 4 - sandy clay loam 5 - clay loam

Default Texture used by DSM and quick and dirty conversions 1 - loamy > loam 2 - fine-loamy > silty clay loam 3 - coarse-loamy > sandy loam

# original before depth class substitution
#reclasssoil<- c(1,11,2,12,3,9,4,5,5,6,6,0,7,0,8,0,9,0,10,0,11,0,12,0,13,0,14,0)

##this has been modified for quick implementation of shallow-deep classes
reclasssoil<- c(1,12,2,26,3,9,4,23,5,5,6,19)


reclasssoilmat<- matrix(reclasssoil, ncol = 2, byrow = TRUE)

reclasssoilmat
##      [,1] [,2]
## [1,]    1   12
## [2,]    2   26
## [3,]    3    9
## [4,]    4   23
## [5,]    5    5
## [6,]    6   19
reclassreproj<- reclassify(masked60, reclasssoilmat)


plot(reclassreproj, main = "SSURGO SDA Top 50cm DC Texture Method", zlim = c(1,14), col = viridis(14))

writeRaster(reclassreproj, file = 'grassexport/cwt_ws32/ssurgosoil.tif', overwrite = TRUE, options=c("COMPRESS=LZW"))

#plot(SoilsMap2, main = "SSURGO USACE Texture Method")



reclassdsm<- c(1,9,2,3,3,12,4,0,5,0,6,0,7,0,8,0,9,0,10,0,11,0,12,0,13,0,14,0)

reclassdsmmat<- matrix(reclassdsm, ncol = 2, byrow = TRUE)

reclassdsmmat
##       [,1] [,2]
##  [1,]    1    9
##  [2,]    2    3
##  [3,]    3   12
##  [4,]    4    0
##  [5,]    5    0
##  [6,]    6    0
##  [7,]    7    0
##  [8,]    8    0
##  [9,]    9    0
## [10,]   10    0
## [11,]   11    0
## [12,]   12    0
## [13,]   13    0
## [14,]   14    0
reclassmaskeddsm <- reclassify(maskeddsm, reclassdsmmat)

plot(reclassmaskeddsm, main = "DSM Top 50cm DC Texture Method", zlim = c(1,14), col = viridis(14))

#plot(maskeddsm)



reclassstatic<- c(1,9,2,9,3,9,4,9,5,9,6,9,7,9,8,9,9,9,10,9,11,9,12,9,13,9,14,9)

reclassstaticmat<- matrix(reclassstatic, ncol = 2, byrow = TRUE)

reclassstaticmat
##       [,1] [,2]
##  [1,]    1    9
##  [2,]    2    9
##  [3,]    3    9
##  [4,]    4    9
##  [5,]    5    9
##  [6,]    6    9
##  [7,]    7    9
##  [8,]    8    9
##  [9,]    9    9
## [10,]   10    9
## [11,]   11    9
## [12,]   12    9
## [13,]   13    9
## [14,]   14    9
reclassmaskedstatic <- reclassify(maskeddsm, reclassstaticmat)

plot(reclassmaskedstatic, main = "Static Loam Texture", zlim = c(1,14), col = viridis(14))

NLCD Land Cover Classification Legend

11 Open Water 12 Perennial Ice/ Snow 21 Developed, Open Space 22 Developed, Low Intensity 23 Developed, Medium Intensity 24 Developed, High Intensity 31 Barren Land (Rock/Sand/Clay) 41 Deciduous Forest 42 Evergreen Forest 43 Mixed Forest 51 Dwarf Scrub 52 Shrub/Scrub 71 Grassland/Herbaceous 72 Sedge/Herbaceous 73 Lichens 74 Moss 81 Pasture/Hay 82 Cultivated Crops 90 Woody Wetlands 95 Emergent Herbaceous Wetlands

Because default RHESSys definition files are being used we will reclassify some of this land cover classification for our landscape

41 Deciduous Forest <- 3 veg_deciduous 42 Evergreen Forest <- 4 veg_evergreen

LANDCOVER<- get_nlcd(template = x.buff.gcs, label = 'coweeta', dataset = "landcover")
#CANOPY<- get_nlcd(template = x.buff.gcs, label = 'coweeta', year = 2016, dataset = "canopy")
IMPERVIOUS<- get_nlcd(template = x.buff.gcs, label = 'coweeta', dataset = "impervious")

plot(LANDCOVER, main = "NLCD Landcover")

#plot(CANOPY)
plot(IMPERVIOUS, main = "NLCD Impervious")

LANDCOVERREPROJ <-  projectRaster(from = LANDCOVER, crs = proj4string(ws32polygon), method = 'ngb')
#CANOPYREPROJ <-  projectRaster(from = CANOPY, crs = proj4string(ws32polygon), method = 'ngb')
IMPERVIOUSREPROJ <-  projectRaster(from = IMPERVIOUS, crs = proj4string(ws32polygon), method = 'ngb')


landcovercrop <- crop(LANDCOVERREPROJ,ws32polygon)
#canopycrop <- crop(CANOPYREPROJ,ws32polygon)
imperviouscrop <- crop(IMPERVIOUSREPROJ,ws32polygon)

maskedlandcover <- mask(landcovercrop,ws32polygon)
#maskedcanopy <- mask(canopycrop,ws32polygon)
maskedimpervious <- mask(imperviouscrop,ws32polygon)


plot(maskedlandcover, main = "Masked Landcover")

#plot(maskedcanopy)
plot(maskedimpervious, main = "Masked Impervious")

{plot(maskedlandcover, main = "Masked Landcover with WS Polygon")
plot(ws32polygon, add = TRUE, col = NA)}

{plot(maskedimpervious, main = "Masked Impervious with WS Polygon")
plot(ws32polygon, add = TRUE, col = NA, border = "white")}

extent(LANDCOVER)
## class      : Extent 
## xmin       : 1128045 
## xmax       : 1136025 
## ymin       : 1402515 
## ymax       : 1411275
extent(LANDCOVERREPROJ)
## class      : Extent 
## xmin       : 271758.1 
## xmax       : 281385.9 
## ymin       : 3875918 
## ymax       : 3886032
extent(maskedlandcover)
## class      : Extent 
## xmin       : 275077 
## xmax       : 276123.5 
## ymin       : 3880946 
## ymax       : 3881563
extent(ws32polygon)
## class      : Extent 
## xmin       : 275066.7 
## xmax       : 276126.7 
## ymin       : 3880940 
## ymax       : 3881570
# template_read(template)
# would be useful to modify or read RHESSys template in R, currently template_read does not display nicely

# Create Templates

{templatetext<- "1
../defs/basin.def
1
../defs/hillslope.def
1
../defs/zone.def
14
../defs/soil_clay.def
../defs/soil_clayloam.def
../defs/soil_loam.def
../defs/soil_loamysand.def
../defs/soil_rock.def
../defs/soil_sand.def
../defs/soil_sandyclay.def
../defs/soil_sandyclayloam.def
../defs/soil_sandyloam.def
../defs/soil_silt.def
../defs/soil_siltyclay.def
../defs/soil_siltyclayloam.def
../defs/soil_siltyloam.def
../defs/soil_water.def
../defs/soil_shallowloam.def
../defs/soil_shallowsandyclayloam
../defs/soil_shallowsandyloam
1
../defs/lu_undev.def
1
../defs/veg_evergreen.def
1
../clim/cwtws32daymet.base
_world basin_ws32.tif 1
 _basin basin_ws32.tif 1
    x       value 0.0
    y       value 0.0
    z       aver dem10.tif
    basin_parm_ID   dvalue 1
    latitude    value 35.04 
    basin_n_basestations    dvalue 0    
    _hillslope subbasins.tif 1
        x       value 0.0
        y       value 0.0
        z       aver dem10.tif
        hill_parm_ID    dvalue 1    
        gw.storage  value 0.0
        gw.NO3      value 0.0
        hillslope_n_basestations    dvalue 0    
        _zone patch.tif 1 
            x         value 0.0
            y         value 0.0
            z         aver dem10.tif 
            zone_parm_ID      dvalue 1  
            area          area  
            slope         aver slope10.tif  
            aspect        spavg aspect10.tif slope10.tif    
            precip_lapse_rate value 1.0
            e_horizon     eqn 0.001 0 e_horizon_1000.tif    
            w_horizon     eqn 0.001 0 w_horizon_1000.tif    
            zone_n_basestations   dvalue 1  
            zone_basestation_ID   dvalue 101    
             _patch patch.tif 1
                x           value 0.0
                y           value 0.0
                z           aver  dem10.tif
                soil_parm_ID            mode ssurgosoil.tif 
                landuse_parm_ID     dvalue 1
                fire_parm_ID        value 1.0
                area            area    
                slope           aver  slope10.tif           
                lna             eqn 0.001 0 topidx10_100.tif    
                Ksat_vertical       value 1.0 
                mpar                value 0.12  
                rz_storage      value 0.0
                unsat_storage       value 0.0
                sat_deficit         value 0.0
                snowpack.water_equivalent_depth value 0.28 
                snowpack.water_depth    value 0.0
                snowpack.T      value -10.0
                snowpack.surface_age    value 0.0
                snowpack.energy_deficit value -0.5 
                litter.cover_fraction   value 1.0
                litter.rain_stored  value 0.00001544
                litter_cs.litr1c    value 0.00001945
                litter_ns.litr1n    value 0.00000082
                litter_cs.litr2c    value 0.00316727
                litter_cs.litr3c    value 0.00021824
                litter_cs.litr4c    value 0.00532178
                soil_cs.soil1c      value 0.00238081
                soil_ns.sminn       value 0.00006014
                soil_ns.nitrate     value 0.00022041
                soil_cs.soil2c      value 0.02335700
                soil_cs.soil3c      value 0.06099644
                soil_cs.soil4c      value 0.37339245
                patch_n_basestations        dvalue 0
                _canopy_strata patch.tif 1  
                    veg_parm_ID     dvalue 3    
                            cover_fraction      value 1.0
                    gap_fraction        value 0.0 
                    rootzone.depth      value 1.0 
                    snow_stored     value 0.0 
                    rain_stored     value 0.0 
                    cs.pool         value 0.0 
                    cs.leafc        value 0.0 
                    cs.dead_leafc       value 0.0 
                    cs.leafc_store      value 0.0 
                    cs.leafc_transfer   value 0.0 
                    cs.live_stemc       value 0.0 
                    cs.livestemc_store  value 0.0 
                    cs.live.stemc_transfer  value 0.0 
                    cs.dead_stem        value 0.0 
                    cs.deadstemc_store  value 0.0 
                    cs.deadstemc_transfer   value 0.0 
                    cs.live_crootc      value 0.0 
                    cs.livecrootc_store value 0.0 
                    cs.livecrootc_transfer  value 0.0 
                    cs.dead_crootc      value 2.0 
                    cs.deadcrootc_store value 0.0 
                    cs.deadcrootc_transfer  value 0.0 
                    cs.frootc       value 0.0 
                    cs.frootc_store     value 0.0 
                    cs.frootc_transfer  value 0.0 
                    cs_cwdc         value 0.0 
                    epv.prev_leafcalloc value 0.0 
                    ns.pool         value 0.0 
                    ns.leafn        value 0.0 
                    ns.dead_leafn       value 0.0 
                    ns.leafn_store      value 0.0 
                    ns.leafn_transfer   value 0.0 
                    ns.live_stemn       value 0.0 
                    ns.livestemn_store  value 0.0 
                    ns.livestemn_transfer   value 0.0 
                    ns.dead_stem        value 0.0 
                    ns.deadstemn_store  value 0.0 
                    ns.deadstemn_transfer   value 0.0 
                    ns.live_crootn      value 0.0 
                    ns.livecrootn_store value 0.0 
                    ns.livecrootn_transfer  value 0.0 
                    ns.dead_crootn      value 0.0 
                    ns.deadcrootn_store value 0.0 
                    ns.deadcrootn_transfer  value 0.0 
                    ns.frootn       value 0.0 
                    ns.frootn_store     value 0.0 
                    ns.frootn_transfer  value 0.0 
                    ns.cwdn         value 0.0 
                    ns.retransn     value 0.0 
                    epv.wstress_days        dvalue 0  
                    epv.min_fparabs     value 0.0 
                    epv.min_vwc     value 0.0 
                    canopy_strata_n_basestations dvalue 0 


"}

write(templatetext,file = "templates/RHESSYSinR_Coweeta_WS32_Spinup.template")

Locate RHESSys Spatial Inputs and Specify Template to be used

Create RHESSys Worldfile and Flowtable for 7.4 based on Spatial Inputs and Template

Raster files and template information are used to setup RHESSys worldfile. Flowtable created by “RHESSysPreprocess” must be converted to work. Unsure why. Preprocessing is controlled by template file.

Template file ‘CWWS32_STATIC.template’ is used to run model normally template file ‘CWWS32_STATICspinup.template’ is used to spin up model. Only difference is artificially extended climate dataset.

#setwd(system.file("extdata/", package = "RHESSysIOinR"))

type = "raster"
typepars = "grassexport/cwt_ws32"
template = "templates/RHESSYSinR_Coweeta_WS32_Spinup.template"
# Suffixes of .world, .flow, and .meta will be appended to the worldfile, flowtable, and metadata files respectively.
name = "CWWS32"
overwrite =  TRUE
streams = "streams.tif"
# Optional Flowtable Spatial Data
roads = "roads.tif"
# impervious = "impervious_map"
# roofs = "roofs_map"
header = TRUE

# With certain configurations roads does not work. Returns message: road width cannot be 0 during preprocessing

RHESSysPreprocess(
  template = template,
  name = name,
  type = type,
  typepars = typepars,
  streams = streams,
  overwrite = overwrite,
  header = header,
  parallel = FALSE,
  make_stream = TRUE)
## [1] Begin world_gen.R
## Reading in maps
## Finished reading in maps
## [1] Writing worldfile
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |======================================================================| 100%
## [1] Created worldfile: ./CWWS32.world
## [1] Created header file: ./CWWS32.hdr
## [1] Begin CreateFlownet.R
## Reading in maps
## Finished reading in maps
## [1] Building flowtable
## [1] Generating unique IDs
## [1] Finding patch neighbors
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |                                                                      |   1%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |==                                                                    |   4%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   6%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=======                                                               |  11%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |=========                                                             |  13%
  |                                                                            
  |=========                                                             |  14%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |==============                                                        |  21%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |================                                                      |  24%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |=================                                                     |  25%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |===================                                                   |  26%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |=====================                                                 |  31%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |======================                                                |  32%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |=======================                                               |  34%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |=========================                                             |  35%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |==========================                                            |  36%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |============================                                          |  39%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |============================                                          |  41%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |=============================                                         |  42%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |==============================                                        |  44%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |===============================                                       |  45%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |=================================                                     |  46%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |==================================                                    |  49%
  |                                                                            
  |===================================                                   |  49%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===================================                                   |  51%
  |                                                                            
  |====================================                                  |  51%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |=====================================                                 |  54%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |========================================                              |  56%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |=========================================                             |  58%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |==========================================                            |  59%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |==========================================                            |  61%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |============================================                          |  62%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |============================================                          |  64%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  66%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |===============================================                       |  68%
  |                                                                            
  |================================================                      |  68%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |=================================================                     |  69%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |=================================================                     |  71%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |===================================================                   |  74%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |=====================================================                 |  75%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |======================================================                |  76%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=======================================================               |  78%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |========================================================              |  79%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |========================================================              |  81%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |=========================================================             |  82%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |==========================================================            |  84%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=============================================================         |  86%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |===============================================================       |  89%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |===============================================================       |  91%
  |                                                                            
  |================================================================      |  91%
  |                                                                            
  |================================================================      |  92%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |=================================================================     |  94%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |====================================================================  |  96%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |====================================================================  |  98%
  |                                                                            
  |===================================================================== |  98%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================|  99%
  |                                                                            
  |======================================================================| 100%
## [1] Buildling flowtable list
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |                                                                      |   1%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |==                                                                    |   4%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   6%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=======                                                               |  11%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |=========                                                             |  13%
  |                                                                            
  |=========                                                             |  14%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |==============                                                        |  21%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |================                                                      |  24%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |=================                                                     |  25%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |===================                                                   |  26%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |=====================                                                 |  31%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |======================                                                |  32%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |=======================                                               |  34%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |=========================                                             |  35%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |==========================                                            |  36%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |============================                                          |  39%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |============================                                          |  41%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |=============================                                         |  42%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |==============================                                        |  44%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |===============================                                       |  45%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |=================================                                     |  46%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |==================================                                    |  49%
  |                                                                            
  |===================================                                   |  49%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===================================                                   |  51%
  |                                                                            
  |====================================                                  |  51%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |=====================================                                 |  54%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |========================================                              |  56%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |=========================================                             |  58%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |==========================================                            |  59%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |==========================================                            |  61%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |============================================                          |  62%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |============================================                          |  64%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  66%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |===============================================                       |  68%
  |                                                                            
  |================================================                      |  68%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |=================================================                     |  69%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |=================================================                     |  71%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |===================================================                   |  74%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |=====================================================                 |  75%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |======================================================                |  76%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=======================================================               |  78%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |========================================================              |  79%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |========================================================              |  81%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |=========================================================             |  82%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |==========================================================            |  84%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=============================================================         |  86%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |===============================================================       |  89%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |===============================================================       |  91%
  |                                                                            
  |================================================================      |  91%
  |                                                                            
  |================================================================      |  92%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |=================================================================     |  94%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |====================================================================  |  96%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |====================================================================  |  98%
  |                                                                            
  |===================================================================== |  98%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================|  99%
  |                                                                            
  |======================================================================| 100%
## [1] Filling 1 pits.
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |======================================================================| 100%
## [1] Writing flowtable
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |                                                                      |   1%
  |                                                                            
  |=                                                                     |   1%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |==                                                                    |   4%
  |                                                                            
  |===                                                                   |   4%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   6%
  |                                                                            
  |=====                                                                 |   7%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |======                                                                |   8%
  |                                                                            
  |======                                                                |   9%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |=======                                                               |  10%
  |                                                                            
  |=======                                                               |  11%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |========                                                              |  12%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |=========                                                             |  13%
  |                                                                            
  |=========                                                             |  14%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |==========                                                            |  15%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |============                                                          |  18%
  |                                                                            
  |=============                                                         |  18%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |==============                                                        |  21%
  |                                                                            
  |===============                                                       |  21%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  22%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |================                                                      |  24%
  |                                                                            
  |=================                                                     |  24%
  |                                                                            
  |=================                                                     |  25%
  |                                                                            
  |==================                                                    |  25%
  |                                                                            
  |==================                                                    |  26%
  |                                                                            
  |===================                                                   |  26%
  |                                                                            
  |===================                                                   |  27%
  |                                                                            
  |===================                                                   |  28%
  |                                                                            
  |====================                                                  |  28%
  |                                                                            
  |====================                                                  |  29%
  |                                                                            
  |=====================                                                 |  29%
  |                                                                            
  |=====================                                                 |  30%
  |                                                                            
  |=====================                                                 |  31%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |======================                                                |  32%
  |                                                                            
  |=======================                                               |  32%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |=======================                                               |  34%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |========================                                              |  35%
  |                                                                            
  |=========================                                             |  35%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |==========================                                            |  36%
  |                                                                            
  |==========================                                            |  37%
  |                                                                            
  |==========================                                            |  38%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |===========================                                           |  39%
  |                                                                            
  |============================                                          |  39%
  |                                                                            
  |============================                                          |  40%
  |                                                                            
  |============================                                          |  41%
  |                                                                            
  |=============================                                         |  41%
  |                                                                            
  |=============================                                         |  42%
  |                                                                            
  |==============================                                        |  42%
  |                                                                            
  |==============================                                        |  43%
  |                                                                            
  |==============================                                        |  44%
  |                                                                            
  |===============================                                       |  44%
  |                                                                            
  |===============================                                       |  45%
  |                                                                            
  |================================                                      |  45%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |=================================                                     |  46%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |=================================                                     |  48%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |==================================                                    |  49%
  |                                                                            
  |===================================                                   |  49%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |===================================                                   |  51%
  |                                                                            
  |====================================                                  |  51%
  |                                                                            
  |====================================                                  |  52%
  |                                                                            
  |=====================================                                 |  52%
  |                                                                            
  |=====================================                                 |  53%
  |                                                                            
  |=====================================                                 |  54%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |======================================                                |  55%
  |                                                                            
  |=======================================                               |  55%
  |                                                                            
  |=======================================                               |  56%
  |                                                                            
  |========================================                              |  56%
  |                                                                            
  |========================================                              |  57%
  |                                                                            
  |========================================                              |  58%
  |                                                                            
  |=========================================                             |  58%
  |                                                                            
  |=========================================                             |  59%
  |                                                                            
  |==========================================                            |  59%
  |                                                                            
  |==========================================                            |  60%
  |                                                                            
  |==========================================                            |  61%
  |                                                                            
  |===========================================                           |  61%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |============================================                          |  62%
  |                                                                            
  |============================================                          |  63%
  |                                                                            
  |============================================                          |  64%
  |                                                                            
  |=============================================                         |  64%
  |                                                                            
  |=============================================                         |  65%
  |                                                                            
  |==============================================                        |  65%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  66%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |===============================================                       |  68%
  |                                                                            
  |================================================                      |  68%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |=================================================                     |  69%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |=================================================                     |  71%
  |                                                                            
  |==================================================                    |  71%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |===================================================                   |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |===================================================                   |  74%
  |                                                                            
  |====================================================                  |  74%
  |                                                                            
  |====================================================                  |  75%
  |                                                                            
  |=====================================================                 |  75%
  |                                                                            
  |=====================================================                 |  76%
  |                                                                            
  |======================================================                |  76%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |======================================================                |  78%
  |                                                                            
  |=======================================================               |  78%
  |                                                                            
  |=======================================================               |  79%
  |                                                                            
  |========================================================              |  79%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |========================================================              |  81%
  |                                                                            
  |=========================================================             |  81%
  |                                                                            
  |=========================================================             |  82%
  |                                                                            
  |==========================================================            |  82%
  |                                                                            
  |==========================================================            |  83%
  |                                                                            
  |==========================================================            |  84%
  |                                                                            
  |===========================================================           |  84%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |============================================================          |  85%
  |                                                                            
  |============================================================          |  86%
  |                                                                            
  |=============================================================         |  86%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |=============================================================         |  88%
  |                                                                            
  |==============================================================        |  88%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |===============================================================       |  89%
  |                                                                            
  |===============================================================       |  90%
  |                                                                            
  |===============================================================       |  91%
  |                                                                            
  |================================================================      |  91%
  |                                                                            
  |================================================================      |  92%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |=================================================================     |  94%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |====================================================================  |  96%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |====================================================================  |  98%
  |                                                                            
  |===================================================================== |  98%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================|  99%
  |                                                                            
  |======================================================================| 100%
## [1] Created flowtable: ./CWWS32.flow
## Make sure to have parallel turned off  when preprocessing, convert flowtables after creation for compatibility
convert_flowtable("CWWS32.flow","CWWS321.flow",make_stream=4)
## [1] "Reading in flow table"
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=======================                                               |  33%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |======================================================================| 100%

Setup RHESSys Inputs/Outputs for Spinup Note that flowtable is converted flowtable. Using original flow table will result in errors.

Original spinup 0219 to 2017. Modified output to 1776-2017 to test code without very long spin-up. Climate files have been modified to original DAYMET timeseries.

## Read climate file to see length available climate data
read_clim("clim/cwtws32extended.rain", dates_out=TRUE)
## [1] "1769 12 27 01" "2021 12 31 24"
read_clim("clim/cwtws32daymet.rain", dates_out=TRUE)
## [1] "1980 01 01 01" "2021 12 31 24"
# Begin RHESSys Spinup
##Short Spinup, original spinup run 1776-2018

input_rhessys = IOin_rhessys_input(
  version = rh_path,
  tec_file = "tecfiles/tec_daily",
  world_file = "CWWS32.world",
  world_hdr_prefix = "CWWS32",
  flowtable = "CWWS321.flow",
  start = "1980 11 1 1",
  end = "2018 11 1 1",
  output_folder = "out",
  output_prefix = "cwws32",
  commandline_options = c("-b -g"))

input_tec_data = IOin_tec_std(start = "1981 11 1 1",
                              end = "2018 11 1 1",
                              output_state = TRUE)

input_hdr = IOin_hdr(
  basin = "defs/basin.def",
  hillslope = "defs/hillslope.def",
  zone = "defs/zone.def",
  soil = c("defs/soil_clay.def","defs/soil_clayloam.def","defs/soil_loam.def","defs/soil_loamysand.def","defs/soil_rock.def","defs/soil_sand.def","defs/soil_sandyclay.def","defs/soil_sandyclayloam.def","defs/soil_sandyloam.def","defs/soil_silt.def","defs/soil_siltyclay.def","defs/soil_siltyclayloam.def","defs/soil_siltyloam.def","defs/soil_water.def", "defs/soil_shallowloam.def", "defs/soil_shallowsandyclayloam.def", "defs/soil_shallowsandyloam.def"),
  landuse = "defs/lu_undev.def",
  stratum = "defs/veg_deciduous.def",
  basestations = "clim/cwtws32daymet.base"
)

##suggested standard hydrologic spinup parameters from rhessys wiki, gw values modified based on previous results from Coweeta Landscape

stdpars<- IOin_std_pars( 
            m = 1,
            k = 10,
            soil_dep= 0.4,
            m_v = 1,
            k_v = 1,
            gw1 = 0.15,
            gw2 = 0.35)

Run Spinup

Spinup runs model on one processor for x number of years. Takes a very long time.

Spinup has been turned off for the purposes of this demonstration. Using previously spun up file ~ 12 hours, run multiple times.

Read Spinup

#spinupoutput<- readin_rhessys_output("out/ws32/cwws32_runspinup")
spinupoutput<- readin_rhessys_output("out/cwws32_runspinup")

spinupoutput$bd$date <- make_date(year=spinupoutput$bd$year, month = spinupoutput$bd$month, day = spinupoutput$bd$day)
spinupoutput$bdg$date <- make_date(year=spinupoutput$bdg$year, month = spinupoutput$bdg$month, day = spinupoutput$bdg$day)

range(spinupoutput$bd$year)
## [1] 1981 2018
plot(spinupoutput$bd$streamflow, type = "l")

plot(spinupoutput$bd$date,spinupoutput$bd$streamflow, type = "l", col = 'black', main = "Predicted Streamflow",
    xlab = "Date", ylab = "Streamflow mm")

plot(spinupoutput$bdg$date, spinupoutput$bdg$lai, type = "l", main = "LAI", xlab = "Date")

{plot(spinupoutput$bdg$soilc~spinupoutput$bdg$date,type = "l", main = "Spinup Soil Carbon", xlab = "Date")
soilctrend <- glm(spinupoutput$bdg$soilc~spinupoutput$bdg$date)
abline(soilctrend, col = 'red')}

{plot(spinupoutput$bdg$date,spinupoutput$bdg$soiln, type = "l", main = "Spinup Soil Nitrogen", xlab = "Date")
soilntrend <- glm(spinupoutput$bdg$soiln~spinupoutput$bdg$date)
abline(soilntrend , col = 'green')}

#prepped<- readin_rhessys_output("G:/Coweeta W32/output/basicspinupws32")

Configure RHESSys Inputs/Outputs for a single run.

Use World Statefile created during spin-up. To run model one time at patch scale.

COMMAND LINE OPTIONS

-b Basin output option. Print out response variables for specified basins. -c Canopy stratum output option. Print out response variables for specified strata. -g Grow option. Try to read in dynamic bgc input data and output dynamic bgc parameters. -h Hillslope output option. Print out response variables for specified hillslopes. -p Patch output option. Print out response variables for specified patches. -r Routing option. Gives name of flow_table to define explicit routing connectivity. Also trigger use of explicit routing over TOPMODEL approach. -c Stratum output option. Print out response variables for specified strata.

read_clim("clim/cwtws32daymet.rain", dates_out = TRUE)
## [1] "1980 01 01 01" "2021 12 31 24"
input_rhessys = IOin_rhessys_input(
  version = rh_path,
  tec_file = "tecfiles/tec_daily",
  world_file = "CWWS32.world.Y2018M10D31H1.state",
  world_hdr_prefix = "CWWS32",
  flowtable = "CWWS321.flow",
  start = "2014 11 1 1",
  end = "2018 11 1 1",
  output_folder = "out",
  output_prefix = "cwws32",
  commandline_options = c("-b -g -p"))
 
# commandline_options = c("-b -g -p 1 128 154914 154914")

## TEC file dictates model output, begin output a year in to allow model SM to stabilize
# do not output_state or worldfile may be overwritten as output is created
input_tec_data = IOin_tec_std(start = "2015 11 1 1",
                              end = "2018 11 1 1",
                              output_state = FALSE)

input_hdr = IOin_hdr(
  basin = "defs/basin.def",
  hillslope = "defs/hillslope.def",
  zone = "defs/zone.def",
  soil = c("defs/soil_clay.def","defs/soil_clayloam.def","defs/soil_loam.def","defs/soil_loamysand.def","defs/soil_rock.def","defs/soil_sand.def","defs/soil_sandyclay.def","defs/soil_sandyclayloam.def","defs/soil_sandyloam.def","defs/soil_silt.def","defs/soil_siltyclay.def","defs/soil_siltyclayloam.def","defs/soil_siltyloam.def","defs/soil_water.def", "defs/soil_shallowloam.def", "defs/soil_shallowsandyclayloam.def", "defs/soil_shallowsandyloam.def"),
  landuse = "defs/lu_undev.def",
  stratum = "defs/veg_deciduous.def",
  basestations = "clim/cwtws32daymet.base")

## Calibrated Values 
#stdpars<- IOin_std_pars(m = 1.12, k = 25.23, soil_dep= 0.35, m_v = 4.07, k_v = 0.59, gw1 = 0.09, gw2 = 0.20)

stdpars<- IOin_std_pars(
            m = 1.66,
            k = 26.16,
            soil_dep= 0.33,
            m_v = 0.31,
            k_v = 11.40,
            gw1 = 0.04,
            gw2 = 0.21)

Run RHESSys once

Read RHESSys Single Run Output

getwd()
## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"
singleruntest<- readin_rhessys_output("out/cwws32_runsingle")

plot(singleruntest$bd$streamflow~singleruntest$bd$date, type = "l", main = "Streamflow", xlab = "Date", ylab = "Streamflow", col = 'DarkBlue')

plot(as.Date(singleruntest$bd$date),singleruntest$bd$rz_storage, type = "l", col = 'black', main = "Basin Scale Root Zone Storage",
    xlab = "Date", ylab = "mm")

#basin scale soil moisture
plot(singleruntest$bd$rz_storage/singleruntest$bdg$root_depth~singleruntest$bd$date, type = "l", main = "RZ_Storage/Root_Depth x Date", xlab = "Date", ylab = "rz_Storage/root_depth", col = 'brown')

plot(singleruntest$bd$lai~singleruntest$bd$date, type = "l", main = "LAI", xlab = "Date", ylab = "LAI", col = 'DarkGreen')

plot(singleruntest$bd$pet~singleruntest$bd$date, type = "l", main = "Potential Evapotranspiration", xlab = "Date", ylab = "PET", col = 'DarkBlue')

plot(singleruntest$bd$et~singleruntest$bd$date, type = "l", main = "Evapotranspiration", xlab = "Date", ylab = "ET", col = 'darkslategray')

plot((singleruntest$bd$unsat_stor/singleruntest$bd$sat_def_z)~singleruntest$bd$date, type = "l", main = "unsat_stor/sat_def_z", xlab = "Date", ylab = "vwc", col = 'DarkGreen')

Plot and Filter Observed Soil Moisture

## Data has not been pre-prepared and must be cleaned

smws32_1<- read.delim("C:/Users/Carlos/Desktop/ORISE/Soil Moisture Data/CW/knb-lter-cwt.1308.19/CWT_132_1_0.txt", header = T)
smws32_2<- read.delim("C:/Users/Carlos/Desktop/ORISE/Soil Moisture Data/CW/knb-lter-cwt.1308.19/CWT_232_1_1308.txt", header = T)
smws32_3<- read.delim("C:/Users/Carlos/Desktop/ORISE/Soil Moisture Data/CW/knb-lter-cwt.1308.19/CWT_332_1_1308.txt", header = T)

smws32_1$Date<- as.Date(smws32_1$Date)
smws32_2$Date<- as.Date(smws32_2$Date)
smws32_3$Date<- as.Date(smws32_3$Date)

smws32_1 <- smws32_1[,!sapply(smws32_1,is.character)]
smws32_2 <- smws32_2[,!sapply(smws32_2,is.character)]
smws32_3 <- smws32_3[,!sapply(smws32_3,is.character)]

#range(smws32_3$smois30A)

smws32_1mean<- aggregate(smws32_1, by = list(smws32_1$Date), FUN = function(x) mean(x, na.rm = TRUE))
smws32_2mean<- aggregate(smws32_2, by = list(smws32_2$Date), FUN = function(x) mean(x, na.rm = TRUE))
smws32_3mean<- aggregate(smws32_3, by = list(smws32_3$Date), FUN = function(x) mean(x, na.rm = TRUE))

smws32_1mean <- smws32_1mean[order(smws32_1mean$Date),]
smws32_2mean <- smws32_2mean[order(smws32_2mean$Date),]
smws32_3mean <- smws32_3mean[order(smws32_3mean$Date),]

#Filter Soil Moisture Datasets

smws32_1mean  <- smws32_1mean[smws32_1mean$smois30A>'0',]
smws32_3mean  <- smws32_3mean[smws32_3mean$smois30A<='0.45',]
smws32_3mean  <- smws32_3mean[smws32_3mean$Date<='2018-05-15',]


plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", main = "Soil Moisture 30cm Site 1", xlab = "Date", ylab = "Soil Moisture 30A")

plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "blue", main = "Soil Moisture 30cm Site 2", xlab = "Date", ylab = "Soil Moisture 30A")

plot(smws32_3mean$Group.1,smws32_3mean$smois30A, type = "l", col = "red", main = "Soil Moisture 30cm Site 3", xlab = "Date", ylab = "Soil Moisture 30A")

{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Observed Soil Moisture WS32 30cm")
lines(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_3mean$Group.1,smws32_3mean$smois30A, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM1","SM2","SM3"), col=c("red","blue","black"), lty=1)}

smws32_1mean  <- smws32_1mean[smws32_1mean$Date>='2017-01-01',]
smws32_2mean  <- smws32_2mean[smws32_2mean$Date>='2017-01-01',]
smws32_3mean  <- smws32_3mean[smws32_3mean$Date>='2017-01-01',]

smws32_1mean  <- smws32_1mean[smws32_1mean$Date<='2018-05-15',]
smws32_2mean  <- smws32_2mean[smws32_2mean$Date<='2018-05-15',]
smws32_3mean  <- smws32_3mean[smws32_3mean$Date<='2018-05-15',]

{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 30cm")
lines(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_3mean$Group.1,smws32_3mean$smois30A, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM1","SM2","SM3"), col=c("red","blue","black"), lty=1)}

{plot(density(smws32_1mean$smois30A), ylim = c(0,14), xlim = c(0,0.35), col = 'red', main = "Density of Observed Soil Moisture WS32 30cm")
lines(density(smws32_2mean$smois30A), col = 'blue')
lines(density(smws32_3mean$smois30A), col = 'black')
legend("topright",inset = 0.05, legend=c("SM1","SM2","SM3"), col=c("red","blue","black"), lty=1)}

{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 Station 1")
lines(smws32_1mean$Group.1,smws32_1mean$smois30B, type = "l", col = "orange", ylim = c(0,0.4))
lines(smws32_1mean$Group.1,smws32_1mean$smois60A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_1mean$Group.1,smws32_1mean$smois60B, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM1-30A","SM1-30B","SM1-60A","SM1-60B"), col=c("red","orange","blue","black"), lty=1)}

{plot(density(smws32_1mean$smois30A), ylim = c(0,14), xlim = c(0,0.35), col = 'red', main = "Density of Observed Soil Moisture WS32 Station 1")
lines(density(smws32_1mean$smois30B), col = 'orange')
lines(density(smws32_1mean$smois60A), col = 'blue')
lines(density(smws32_1mean$smois60B), col = 'black')
legend("topright",inset = 0.05, legend=c("SM1-30A","SM1-30B","SM1-60A","SM1-60B"), col=c("red","orange","blue","black"), lty=1)}

{plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 Station 2")
lines(smws32_2mean$Group.1,smws32_2mean$smois30B, type = "l", col = "orange", ylim = c(0,0.4))
lines(smws32_2mean$Group.1,smws32_2mean$smois60A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_2mean$Group.1,smws32_2mean$smois60B, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM2-30A","SM2-30B","SM2-60A","SM2-60B"), col=c("red","orange","blue","black"), lty=1)}

##extra filtering for smws32_60a
smws32_2mean <- smws32_2mean %>%
  mutate(smois60A = replace(smois60A,Date >= "2017-11-26", NA ))

{plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 Station 2 Filtered")
lines(smws32_2mean$Group.1,smws32_2mean$smois30B, type = "l", col = "orange", ylim = c(0,0.4))
lines(smws32_2mean$Group.1,smws32_2mean$smois60A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_2mean$Group.1,smws32_2mean$smois60B, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM2-30A","SM2-30B","SM2-60A","SM2-60B"), col=c("red","orange","blue","black"), lty=1)}

{plot(smws32_3mean$Group.1,smws32_3mean$smois30A, type = "l", col = "red", ylim = c(0,0.4), xlab = "Date", ylab = "Soil Moisture", main = "Truncated Observed Soil Moisture WS32 Station 3")
lines(smws32_3mean$Group.1,smws32_3mean$smois30B, type = "l", col = "orange", ylim = c(0,0.4))
lines(smws32_3mean$Group.1,smws32_3mean$smois60A, type = "l", col = "blue", ylim = c(0,0.4))
lines(smws32_3mean$Group.1,smws32_3mean$smois60B, type = "l", col = "black", ylim = c(0,0.4))
legend("topright",inset = 0.05, legend=c("SM3-30A","SM3-30B","SM3-60A","SM3-60B"), col=c("red","orange","blue","black"), lty=1)}

soilmoisture30a<- data.frame(smws32_1mean$Date,smws32_1mean$smois30A, smws32_1mean$smois30B,smws32_1mean$smois60A,smws32_1mean$smois60B)
soilmoisture30b<- data.frame(smws32_2mean$Date,smws32_2mean$smois30A, smws32_2mean$smois30B,smws32_2mean$smois60A,smws32_2mean$smois60B)
soilmoisture30c<- data.frame(smws32_3mean$Date,smws32_3mean$smois30A,smws32_3mean$smois30B,smws32_3mean$smois60A,smws32_3mean$smois60B)


merged30a<- merge(soilmoisture30a,soilmoisture30b, by.x = "smws32_1mean.Date", by.y = "smws32_2mean.Date", all = FALSE)

merged30b<- merge(merged30a,soilmoisture30c, by.x = "smws32_1mean.Date", by.y = "smws32_3mean.Date", all = TRUE)

merged30c<- as.data.frame(rowMeans((merged30b[-1]), na.rm= TRUE))


mergedsm<- cbind(merged30b,merged30c)

colnames(mergedsm)<- c("Date","smois30site1a","smois30site1b","smois60site1a","smois60site1b","smois30site2a","smois30site2b","smois60site2a","smois60site2b","smois30site3a","smois30site3b","smois60site3a","smois60site3b","mergedsoilmoisture")

plot(mergedsm$Date, mergedsm$mergedsoilmoisture, type = "l", main = "Average WS32 Soil Moisture", xlab = "Date", ylab = "Merged Soil Moisture All Data")

plot(mergedsm)

## Calculate Calibration/Validation Split based on SM data

# range(mergedsm$Date)

daysofsmdata<- as.numeric(range(mergedsm$Date)[2] - range(mergedsm$Date)[1])

#70/30 split

Caldays<- round(0.7*daysofsmdata)
Valdays<- daysofsmdata-Caldays

Caldates<- as.Date(1)

Caldates[1]<- range(mergedsm$Date)[1]

Caldates[2]<- range(mergedsm$Date)[1]+Caldays

Valdates<- as.Date(1)

Valdates[1]<-  range(Caldates)[2]+1

Valdates[2]<- range(mergedsm$Date)[2]

Plot and Filter Observed Streamflow Data

Original Datasets are recorded as Water Year, We convert to calendar year

## import daily flow dataset
obsws32<- read.csv('C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32.csv', header = F) 

# dataset is recorded as water year, here we convert from water year and note that the water year for this landscape begins in November
## calendar year seems to be correct

#colnames(obsws32)<- c("WYR","Month","Day","flow")

obsws32$WYR <- as.numeric(obsws32$V1)
obsws32$Month <- as.numeric(obsws32$V2)
obsws32$Day <- as.numeric(obsws32$V3)
obsws32$flow <- as.numeric(obsws32$V4)

obsws32$CalendarYear<- ifelse((obsws32$Month=="12"|obsws32$Month=="11"),obsws32$WYR-1,obsws32$WYR)

obsws32$CalendarDate <- as.Date(with(obsws32,paste(CalendarYear,Month,Day,sep="-")),"%Y-%m-%d") 

##this removal of 0s may be reason for lag, 0s may be purposefully included in dataset
# instead create timeseries of date from start to end of timeseries and then merge
#obsws32<- obsws32[!(obsws32$flow=="0"),]

obsws32 <- drop_na(obsws32)



plot(obsws32$CalendarDate, obsws32$flow, type = "l", ylab = "Observed Discharge mm/d", xlab = "Date", main = "Observed Discharge WS32 Older Dataset")

obsws32_1<- read.csv('C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2017.csv', header = T)
obsws32_2<- read.csv('C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2018.csv', header = T)
obsws32_3<- read.csv('C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2019.csv', header = T)
obsws32_4<- read.csv('C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2020.csv', header = T)
obsws32_5<- read.csv('C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/daily_flow_ws32_2021.csv', header = T)

obsws32b<- rbind(obsws32_1,obsws32_2,obsws32_3,obsws32_4,obsws32_5)

obsws32b$CalendarDate <- as.Date(with(obsws32b,paste(year,month,day,sep="-")),"%Y-%m-%d")

plot(obsws32b$CalendarDate,obsws32b$flow, type = "l", xlab = "Date", ylab = "Streamflow", main = "Observed Discharge WS32 Newer Datasets Merged")

dataset1   <- data.frame(Date=obsws32$CalendarDate,flow=obsws32$flow)
dataset2   <- data.frame(Date=obsws32b$CalendarDate,flow=obsws32b$flow)

fullset<- rbind(dataset1,dataset2)

### Streamflow at Coweeta logged as CSM
### conversion from EcohydRology project
## because we are converting from CSM we set watershed area to 1

fullset$discharge_mm<- ConvertFlowUnits(cfs = fullset$flow, WA=1, AREAunits = "mi2")

#{plot(fullset$Date,fullset$flow, type = "l", xlab = "Date", ylab = "Flow", main = "Streamflow WS32")
#lines(fullset$Date,fullset$discharge_mm, col = 'red')}

plot(fullset$Date,fullset$discharge_mm, type = "l", xlab = "Date", ylab = "Discharge mm/d", main = "Observed Discharge WS32 All Datasets Merged")

plot(fullset$Date,fullset$discharge_mm, type = "l", xlab = "Date", ylab = " Log Discharge mm/d", main = "Log of Observed Discharge WS32 All Datasets Merged", log = "y")

fullset$observedstreamflow<- fullset$discharge_mm

write.csv(fullset,"C:/Users/Carlos/Desktop/ORISE/Streamflow Data/CW/Streamflow/Daily Streamflow Watershed 32/Daily Streamflow Watershed 32/ObservedWS32.csv")



obsws32<- fullset

obsws32sm<- merge(obsws32,mergedsm, by = "Date")

#Caldates
#Valdates

obsws32cal <- obsws32[obsws32$Date >= Caldates[1] & obsws32$Date <= Caldates[2], ]
  
obsws32val <- obsws32[obsws32$Date >= Valdates[1] & obsws32$Date <= Valdates[2], ]
  
obsws32smcal <- obsws32sm[obsws32sm$Date >= Caldates[1] & obsws32sm$Date <= Caldates[2], ]

obsws32smval <- obsws32sm[obsws32sm$Date >= Valdates[1] & obsws32sm$Date <= Valdates[2], ]

Merge Observed and Simulated Data

## merge values instead of subsetting to match up with simulated dates
singleruntest$bd<- merge(singleruntest$bd,obsws32, by.x = "date", by.y = "Date", all = FALSE)

## merge soil moisture values but include all data
singleruntest$data <-merge(singleruntest$bd,obsws32sm, by.x = "date", by.y = "Date", all = TRUE)

Plot Observed and Simulated Data, Plot RHESSys Outputs

## Plot Hydrograph comparing modeled and observed streamflow
{hydrograph(input=singleruntest$bd, streamflow=singleruntest$bd$observedstreamflow, streamflow2 = singleruntest$bd$streamflow, timeSeries = singleruntest$bd$date, precip = singleruntest$bd$streamflow,
           P.units = "mm", S.units = "mm normalized by basin area", S1.col = 'Blue', S2.col = 'Red')

legend("topleft",inset = .2, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("blue","red"), lty=1:2)}

{plot(singleruntest$data$date,singleruntest$data$observedstreamflow.x, type = "l", xlab = "Date", ylab = "Streamflow", main = "Predicted vs Observed Streamflow")
lines(singleruntest$data$date, singleruntest$data$streamflow, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:1)}

{plot(singleruntest$data$date,singleruntest$data$observedstreamflow.x, type = "l", xlab = "Date", ylab = "Streamflow", main = "Predicted vs Observed Log Streamflow", log = "y")
lines(singleruntest$data$date, singleruntest$data$streamflow, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:1)}

predictedvsobservedstreamlm<- lm(observedstreamflow~streamflow, data = singleruntest$bd)

{plot(singleruntest$bd$observedstreamflow,singleruntest$bd$streamflow, xlab = "Observed", ylab = "Predicted", main = "Predicted vs Observed Streamflow", xlim = c(0,55), ylim = c(0,55))
abline(a=0, b=1, col = 'RED', lty = 2, lwd = 2)
legend("right",inset = 0.01, legend = paste("R2 =",format(summary(predictedvsobservedstreamlm)$r.squared,digits=3)))}

{plot(singleruntest$data$date,(singleruntest$data$rz_storage/singleruntest$data$rootdepth), type = "l", ylim = c(0,0.35), xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Soil Moisture")
lines(singleruntest$data$date, singleruntest$data$mergedsoilmoisture, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Soil Moisture","Predicted Soil Moisture"), col=c("black","red"), lty=1:1)}

predictedvsobservedsoillm<- lm((rz_storage/rootdepth)~mergedsoilmoisture, data = singleruntest$data)

{plot((singleruntest$data$rz_storage/singleruntest$data$rootdepth),singleruntest$data$mergedsoilmoisture, xlab = "Observed", ylab = "Predicted", main = "Predicted vs Observed Basin Scale Soil Moisture", xlim = c(0.05,0.4), ylim = c(0.05,0.4))
abline(a=0, b=1, col = 'RED', lty = 2, lwd = 2)
legend("right",inset = 0.01, legend = paste("R2 =",format(summary(predictedvsobservedsoillm)$r.squared,digits=3)))}

## Plot other model outputs

par(mfrow=c(5,1), mar = c(2,5,4,2))

plot(singleruntest$bd$date,singleruntest$bd$tmin, type = "l", col = 'BLUE', ylim = c(-15,30), ylab = "Temperature", xlab = "Date")
lines(singleruntest$bd$date,singleruntest$bd$tmax, type = "l", col = 'RED')
plot(singleruntest$bd$date,singleruntest$bd$precip, type = "h", ylab = "Precipitation", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$streamflow, type = "l", ylab = "Streamflow", xlab = "Date", col = 'blue')
plot(singleruntest$bd$date,singleruntest$bd$baseflow, type = "l", ylab = "Baseflow", xlab = "Date", col = 'brown')
plot(singleruntest$bd$date,singleruntest$bd$rz_storage, type = "l", ylab = "RZ Storage", xlab = "Date", col = 'dark green')

plot(singleruntest$bd$date,singleruntest$bd$psn ,type = "l", ylab = "PSN", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$et, type = "l", ylab = "Evaoptranspiration", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$pet, type = "l", ylab = "PET", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$gw.Qout, type = "l", ylab = "Q out", xlab = "Date")
plot(singleruntest$bd$date,singleruntest$bd$gw.storage, type = "l", ylab = "GW Storage", xlab = "Date")

par(mfrow=c(1,1))

### will only plot in growth mode , used to check that vegetation is initialized properly ##
plot(singleruntest$bdg$date, singleruntest$bdg$lai, type = "l", main = "LAI", xlab = "Date")

plot(singleruntest$bdg$date, singleruntest$bdg$soilc, type = "l", main = "Soil Carbon", xlab = "Date")

## plot cal/val sets

{plot(singleruntest$data$date,(singleruntest$data$rz_storage/singleruntest$data$rootdepth), type = "l", ylim = c(0,0.35), xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Soil Moisture")
lines(obsws32smcal$Date,obsws32smcal$mergedsoilmoisture, col = "BLUE")
lines(obsws32smval$Date,obsws32smval$mergedsoilmoisture, col = "RED")
legend("topright",inset = 0, legend=c("Predicted Soil Moisture","Observed Calibration Soil Moisture", "Observed Validation Soil Moisture"), col=c("black","blue","red"), lty=c(1,1,1))}

{plot(singleruntest$data$date,singleruntest$data$streamflow, type = "l", xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Streamflow")
lines(obsws32cal$Date,obsws32cal$observedstreamflow, col = "BLUE")
lines(obsws32val$Date,obsws32val$observedstreamflow, col = "RED")
legend("topright",inset = 0, legend=c("Predicted Streamflow","Observed Calibration Streamflow", "Observed Validation Streamflow"), col=c("black","blue","red"), lty=c(1,1,1))}

Testing Precip input vs model precip

Spatial Data Display

setwd(system.file("extdata/", package = "RHESSysIOinR"))

# Original input maps have been cropped before processing in R, without cropping the RHESSYs preparation in the R environment will not function properly, below are the input maps used by R
ws32<- raster("grassexport/cwt_ws32/basin_ws32.tif")
patches<- raster("grassexport/cwt_ws32/patch.tif")
acc10<- raster("grassexport/cwt_ws32/acc10.tif")
aspect10<- raster("grassexport/cwt_ws32/aspect10.tif")
#bt1000<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/b.t1000.tif")
dem10<- raster("grassexport/cwt_ws32/dem10.tif")
dem10f<- raster("grassexport/cwt_ws32/dem10f.tif")
dir10<- raster("grassexport/cwt_ws32/dir10.tif")
drain10<- raster("grassexport/cwt_ws32/drain10.tif")
#ht1000<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/h.t1000.tif")
hillslope<- raster("grassexport/cwt_ws32/hillslope.tif")
roads<- brick("grassexport/cwt_ws32/roads.tif")
#shade10<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/shade10.tif")
slope10<- raster("grassexport/cwt_ws32/slope10.tif")
streams<- raster("grassexport/cwt_ws32/streams.tif")
topidx10<- raster("grassexport/cwt_ws32/topidx10_100.tif")
#xmap<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/xmap.tif")
#ymap<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/ymap.tif")


plot(ws32, main ="ws32")

plot(patches, main = "patches")

plot(acc10, main = "acc10")

plot(aspect10, main = "aspect10")

#levelplot(bt1000, col.regions = brewer.pal(name='Spectral', n = 11),at=seq(0,6500,1))
#plot(bt1000, main = "bt1000")
plot(dem10, main = "dem10")

plot(dem10f, main = "dem10f")

plot(dir10, main = "dir10")

plot(drain10, main = "drain10")

#plot(ht1000, main = "ht1000")
plot(hillslope, main = "hillslope")

plot(roads, main = "roads")

#plot(shade10, main = "shade10")
plot(slope10, main = "slope10")

plot(streams, main = "streams")

plot(topidx10, main = "topidx10")

#plot(xmap, main = "xmap")
#plot(ymap, main = "ymap")
ObservedSMLocationsnt <- read_sf("C:/Users/Carlos/Desktop/ORISE/Soil Moisture Data/CW/knb-lter-cwt.1308.19/1308.kml")

ObservedSMLocations <- st_transform(ObservedSMLocationsnt, crs(dem10))

colorramp<- viridis(128)

{plot(dem10, col = colorramp, main = "Coweeta Observed Soil Moisture Locations")
points(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1]), c(ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", lwd = 1)}

ObservedSMLocations$Name[3]
## [1] "132 plot center"
ObservedSMLocations[[3]][[3]]
## POINT (275432.6 3881522)
ObservedSMLocations$Name[4]
## [1] "232 plot center"
ObservedSMLocations[[3]][[4]]
## POINT (275578.2 3881429)
ObservedSMLocations$Name[5]
## [1] "332 plot center"
ObservedSMLocations[[3]][[5]]
## POINT (275760.9 3881284)
matrix(c(ObservedSMLocations$Name[3],ObservedSMLocations$Name[4],ObservedSMLocations$Name[5],ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), nrow = 3, ncol = 3)
##      [,1]              [,2]               [,3]              
## [1,] "132 plot center" "275432.649744231" "3881522.05364127"
## [2,] "232 plot center" "275578.214010129" "3881428.55011708"
## [3,] "332 plot center" "275760.864807773" "3881284.18476691"
## Extract values for single date
displayday <- singleruntest$pd[which(singleruntest$pd$date=="2018-03-21")]

## reclassify all patches to na and then reclassify to model output values
reclass_all<- c(-2147483648, 2147483647,NA)
reclass_allm<- matrix(reclass_all, ncol = 3, byrow = TRUE)
displaydaypatches<- reclassify(patches,reclass_allm)

displayday$calculatedrzsm <- displayday$rz_storage/displayday$root.depth

reclass_displayday <- cbind(displayday$patchID,displayday$calculatedrzsm)
reclass_displayday <- as.matrix(reclass_displayday)

displaydaypatches<- reclassify(patches,reclass_displayday)

## Plot Spatial Data Output for RZ Storage
#plot(displaydaypatches, xlim = c(280900,282500), ylim = c(4870000,4872000))

{plot(mask(displaydaypatches,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 RZ Soil Moisture Prediction for 1 day", line = -2)
points(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1]), c(ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", cex = 2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}

reclass_displayday <- cbind(displayday$patchID,displayday$root.depth)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches2<- reclassify(patches,reclass_displayday)

{plot(mask(displaydaypatches2,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 root.depth Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}

reclass_displayday <- cbind(displayday$patchID,displayday$rz_storage)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches3<- reclassify(patches,reclass_displayday)

{plot(mask(displaydaypatches3,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 rz_storage Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}

reclass_displayday <- cbind(displayday$patchID,displayday$sat_def)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches4<- reclassify(patches,reclass_displayday)

{plot(mask(displaydaypatches4,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 saturation deficit Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}

## Plot Hydrograph
{hydrograph(input=singleruntest$bd, streamflow=singleruntest$bd$observedstreamflow, streamflow2 = singleruntest$bd$streamflow, timeSeries = singleruntest$bd$date, precip = singleruntest$bd$streamflow,
           P.units = "mm", S.units = "mm normalized by basin area", S1.col = 'Blue', S2.col = 'Red')
legend("topleft",inset = .2, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("blue","red"), lty=1:2)}

### end show patch map




patchnumbers<- raster::extract(patches,matrix(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), nrow = 3, ncol = 2))


ObserveSoilMoisturePatches <- matrix(c(ObservedSMLocations$Name[3],ObservedSMLocations$Name[4],ObservedSMLocations$Name[5],ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2],patchnumbers), nrow = 3, ncol = 4)


## Extract values for a single patch, select patch ID
ObserveSoilMoisturePatches[1,4]
## [1] "136557"
displaypatch<- singleruntest$pd[which(singleruntest$pd$patchID==ObserveSoilMoisturePatches[1,4])]
ObserveSoilMoisturePatches[2,4]
## [1] "141942"
displaypatch2<- singleruntest$pd[which(singleruntest$pd$patchID==ObserveSoilMoisturePatches[2,4])]
ObserveSoilMoisturePatches[3,4]
## [1] "149478"
displaypatch3<- singleruntest$pd[which(singleruntest$pd$patchID==ObserveSoilMoisturePatches[3,4])]

plot(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", main = 'CW WS27 More Detailed RZ Storage Patch #147322', ylim = c(0,.4))

{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", ylim = c(0,.4), col = 'darkorange', xlim=as.Date(c("2016-10-1","2019-10-1")), main = "CW WS32 More Detailed Soil Moisture vs Patches",
     xlab = "Date", ylab = "Soil Moisture mm")
lines(smws32_2mean$Group.1,smws32_2mean$smois30A,type = "l", col = 'darkred')
lines(smws32_3mean$Group.1,smws32_3mean$smois30A,type = "l", col = 'blue')
lines(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'red', lty = 2)
lines(displaypatch2$date,(displaypatch2$rz_storage/displaypatch2$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'green', lty = 2)
lines(displaypatch3$date,(displaypatch3$rz_storage/displaypatch3$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'black', lty = 2)

legend("bottomleft",inset = .01, legend=c("Observed Patch 1","Predicted Patch 1","Observed Patch 2","Predicted Patch 2","Observed Patch 3","Predicted Patch 3"), col=c("darkorange", "red","darkred","green","blue","black"), lty=c(1,2,1,2,1,2))}

{plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", ylim = c(0,0.4), col = 'BLUE', xlim=as.Date(c("2016-10-1","2019-10-1")), main = "CW WS32 More Detailed Soil Moisture vs Patch #154914",
     xlab = "Date", ylab = "Soil Moisture mm")
lines(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", main = 'CW WS27 More Detailed RZ Storage Patch #154914', col = 'red')
legend("bottomleft",inset = .01, legend=c("Observed Soil Moisture","Predicted Soil Moisture"), col=c("blue","red"), lty=1:1)}

#
#obsws32cal <- obsws32[obsws32$Date >= Caldates[1] & obsws32$Date <= Caldates[2], ]
#obsws32val <- obsws32[obsws32$Date >= Valdates[1] & obsws32$Date <= Valdates[2], ]
#obsws32smcal <- obsws32sm[obsws32sm$Date >= Caldates[1] & obsws32sm$Date <= Caldates[2], ]
#obsws32smval <- obsws32sm[obsws32sm$Date >= Valdates[1] & obsws32sm$Date <= Valdates[2], ]

head(obsws32cal)
##             Date flow discharge_mm observedstreamflow
## 22961 2017-01-01 2.30     2.172645           2.172645
## 22962 2017-01-02 3.04     2.871669           2.871669
## 22963 2017-01-03 2.72     2.569388           2.569388
## 22964 2017-01-04 2.11     1.993165           1.993165
## 22965 2017-01-05 1.92     1.813686           1.813686
## 22966 2017-01-06 1.81     1.709777           1.709777
head(obsws32val)
##             Date flow discharge_mm observedstreamflow
## 23311 2017-12-17 2.14     2.021504           2.021504
## 23312 2017-12-18 2.27     2.144306           2.144306
## 23313 2017-12-19 2.23     2.106521           2.106521
## 23314 2017-12-20 3.39     3.202289           3.202289
## 23315 2017-12-21 2.74     2.588281           2.588281
## 23316 2017-12-22 2.60     2.456033           2.456033
head(obsws32smcal)
##         Date flow discharge_mm observedstreamflow smois30site1a smois30site1b
## 1 2017-01-01 2.30     2.172645           2.172645     0.2194375     0.2468125
## 2 2017-01-02 3.04     2.871669           2.871669     0.2241528     0.2509896
## 3 2017-01-03 2.72     2.569388           2.569388     0.2184132     0.2409236
## 4 2017-01-04 2.11     1.993165           1.993165     0.1952743     0.2236250
## 5 2017-01-05 1.92     1.813686           1.813686     0.1825451     0.2149410
## 6 2017-01-06 1.81     1.709777           1.709777     0.1745729     0.2097500
##   smois60site1a smois60site1b smois30site2a smois30site2b smois60site2a
## 1     0.1988021     0.2463715     0.2162465     0.2141181     0.2191528
## 2     0.2243646     0.2737917     0.2244410     0.2392847     0.2239028
## 3     0.2177813     0.2657361     0.2209722     0.2291771     0.2154028
## 4     0.1993090     0.2483194     0.2058437     0.2039549     0.1967431
## 5     0.1845556     0.2335799     0.1972882     0.1883021     0.1870486
## 6     0.1744132     0.2231181     0.1914861     0.1774236     0.1808646
##   smois60site2b smois30site3a smois30site3b smois60site3a smois60site3b
## 1     0.1775868     0.2830729     0.2573576     0.2625312     0.2399167
## 2     0.2186528     0.2833924     0.2626007     0.2942917     0.2684410
## 3     0.2180486     0.2697604     0.2480833     0.2926944     0.2704132
## 4     0.1968403     0.2399167     0.2206285     0.2794965     0.2609514
## 5     0.1786215     0.2264861     0.2081111     0.2668299     0.2512847
## 6     0.1669340     0.2187014     0.2009896     0.2569028     0.2437292
##   mergedsoilmoisture
## 1          0.2317839
## 2          0.2490255
## 3          0.2422839
## 4          0.2225752
## 5          0.2099661
## 6          0.2015738
head(obsws32smval)
##           Date flow discharge_mm observedstreamflow smois30site1a smois30site1b
## 351 2017-12-17 2.14     2.021504           2.021504     0.1181493     0.2395313
## 352 2017-12-18 2.27     2.144306           2.144306     0.1165625     0.2392708
## 353 2017-12-19 2.23     2.106521           2.106521     0.1135139     0.2332326
## 354 2017-12-20 3.39     3.202289           3.202289     0.1484514     0.2609861
## 355 2017-12-21 2.74     2.588281           2.588281     0.1391528     0.2532083
## 356 2017-12-22 2.60     2.456033           2.456033     0.1312049     0.2428194
##     smois60site1a smois60site1b smois30site2a smois30site2b smois60site2a
## 351     0.1677882     0.2054826     0.1575382     0.1352778            NA
## 352     0.1635625     0.2064132     0.1605903     0.1370000            NA
## 353     0.1600174     0.2057812     0.1585660     0.1363750            NA
## 354     0.1844479     0.2319132     0.1872778     0.1721458            NA
## 355     0.2092639     0.2517917     0.1925139     0.1736111            NA
## 356     0.1954965     0.2400937     0.1847604     0.1649479            NA
##     smois60site2b smois30site3a smois30site3b smois60site3a smois60site3b
## 351    0.09663194     0.2347188     0.1869757     0.2350000     0.2169896
## 352    0.09784375     0.2385694     0.1832847     0.2350000     0.2167847
## 353    0.09828472     0.2289757     0.1792153     0.2350000     0.2160000
## 354    0.12405208     0.2609826     0.2203681     0.2421354     0.2173993
## 355    0.14059028     0.2525278     0.2136979     0.2712326     0.2365937
## 356    0.13011806     0.2397049     0.1995382     0.2683160     0.2404201
##     mergedsoilmoisture
## 351          0.1812803
## 352          0.1813529
## 353          0.1786329
## 354          0.2045600
## 355          0.2121985
## 356          0.2034018
#subset Validation range based on Val Dates set by available SM data

validationsubset <- singleruntest$bd[singleruntest$bd$date >= Valdates[1] & singleruntest$bd$date <= Valdates[2], ]

## Fit tests for streamflow

NSE(validationsubset$streamflow,validationsubset$observedstreamflow)
## [1] 0.5994421
NSE(validationsubset$streamflow,validationsubset$observedstreamflow, FUN = log, epsilon = "Pushpalatha2012", na.rm=TRUE)
## [1] 0.722282
KGE(validationsubset$streamflow,validationsubset$observedstreamflow)
## [1] 0.802176
plot(validationsubset$streamflow~validationsubset$date, type = "l", col = 'red')
lines(validationsubset$observedstreamflow~validationsubset$date)

# Plot Observed and Predicted Streamflow

{plot(validationsubset$date,validationsubset$streamflow, type = 'l', lty = 2, col = 'RED', main = "Streamflow")
lines(validationsubset$date,validationsubset$observedstreamflow, type = 'l', col = 'BLACK')
legend("topright",inset = .1, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:2)}

{plot(singleruntest$bd$date,singleruntest$bd$streamflow, type = 'l', lty = 2, col = 'RED', main = "Streamflow")
lines(singleruntest$bd$date,singleruntest$bd$observedstreamflow, type = 'l', col = 'BLACK')
legend("topright",inset = .1, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:2)}

#Plot Model RZ SM outputs for patch
{plot(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l",ylim = c(0,0.50), ylab = "Root Zone Storage (mm/m) ", xlab = "Date", main = 'Root Zone Display Patch')
lines(displaypatch$date,(displaypatch$rz_field_capacity/displaypatch$root.depth), type = "l", col = "blue")
lines(displaypatch$date,(displaypatch$rz_wilting_point/displaypatch$root.depth), type = "l", col = "red")
legend("bottomleft",inset = .1, legend=c("Predicted RZ Storage","RZ Field Capacity","RZ Wilting Point"), col=c("black","blue","red"), lty=1:1)}

Setup RHESSys Calibration, CLHS used to reduce amount of samples processed

The model calibrates over the following parameters.

M - Decay of hydraulic conductivity with depth K - Hydraulic conductivity at the surface Sd - soil depth Mv - Vertical decay of hydraulic conductivity with depth Kv - Vertical hydraulic conductivity at the surface GW1 - Saturated to Groundwater Coefficient GW2 - Groundwater Loss Coefficient

Parameter & Typical Range from RHESSys Wiki

M - .01 - 20 K - 1 - 600 GW1 - .001 - .3 GW2 - .01 - .9

" In order to generate an adequate sample across the full range for each parameter, values from .001 - .01, values from .01 - 1, and values over 1 should be generated separately. The m and K parameters are multipliers on the initial values set for m and K (values assigned to the m and K maps)." - Taken from RHESSys Wiki

setwd(system.file("extdata/", package = "RHESSysIOinR"))
# number of samples for full data space
fs <- 10000
# number of samples to calibrate over, 100 samples take approx 40 mins in base R and 1 hour in rmarkdown
n <- 150
# create data frame
d <- data.frame(
  m = runif(fs, min=0.5, max=3),
  k = runif(fs, min=0.5, max=60),
  soil_dep = runif(fs, min=0.1, max=0.5),
  m_v = runif(fs, min=0.5, max=15),
  k_v = runif(fs, min=0.5, max=15),
  gw1 = runif(fs, min=0.1, max=0.3),
  gw2 = runif(fs, min=0.1, max=0.5))
# marginal distributions, for analysis
m <- melt(d)
## No id variables; using all as measure variables
bwplot(variable ~ value, data = m, par.settings = tactile.theme(), main = "distribution of calibration parameters before clhs subsetting")

# cLHS subset
s <- clhs(d, size = n, simple = TRUE, use.cpp = TRUE)
# combine full dataset + subset
g <- make.groups(full = d, cLHS = d[s, ])
# create new list
params<- d[s, ]

params[1:5,]
##              m        k  soil_dep       m_v        k_v       gw1       gw2
## 1152 0.5114497 21.44952 0.1439679 9.6779017  9.5969864 0.1829364 0.4490044
## 1798 0.9612092 16.05240 0.4167406 4.8633875 14.2459499 0.2539121 0.2319639
## 5984 1.3464745 52.34296 0.4217323 1.3377789 14.4170302 0.2133951 0.4848750
## 7019 0.5357427 57.27382 0.1720589 0.5526789  6.8889040 0.1237745 0.3778741
## 8216 2.8337005 51.29474 0.1356276 2.2444764  0.6041038 0.1726981 0.2726171
## Quick and dirty way to include "uncalibrated" result to test against, groundwater model still included for now, unsure of what to do with parameters for gw
## params[1,] = c(1,1,1,1,1,0,0)
## did not produce significantly different result

#write param table to specific text file if file does not exist
## if name is not changed then old table is used
# paramtablename = "out/cwparamtable021723.txt" 

####WORKING

### WORKING

### WORKING
getwd()
## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"
write.table(params, file = "out/cwparamtable032723.txt")

### writes new params table with todays date to a text file
#write.table(params, file = paste("out/cwparamtable",format(Sys.Date(),"%m%d%y"),".txt",sep = ""))

## reads in specific params file, used for keeping input params constant for analysis
## do not rewrite file if it exits

#{if (file.exists("out/cwparamtable102522.txt"))
#   print("text")}


 #params <- read.table(file = "out/cwparamtable102522.txt")
params <- read.table(file = "out/cwparamtable032823.txt")


# prepare RHESSys for run
params<- as.data.frame(params)

input_rhessys = IOin_rhessys_input(
  version = rh_path,
  tec_file = "tecfiles/tec_daily",
  world_file = "CWWS32.world.Y2018M10D31H1.state",
  world_hdr_prefix = "CWWS32",
  flowtable = "CWWS321.flow",
  start = "2014 11 1 1",
  end = "2018 11 1 1",
  output_folder = "out",
  output_prefix = "cwws32",
  commandline_options = c("-b -g "))

#commandline_options = c("-b -p 1 128 154914 154914 ")) <- Patch Specific Output
  
## must incorporate single patch output or local patch output for the purposes of calibration
## make sure to run calibration with no modifiers or 1 for everything except gw1 and gw2
## must refine calibration to do multiple steps, first gw1 gw2, then other variables
## multi point calibration and validation

input_tec_data = IOin_tec_std(start = "2015 11 1 1",
                              end = "2018 11 1 1",
                              output_state = FALSE)

input_hdr = IOin_hdr(
  basin = "defs/basin.def",
  hillslope = "defs/hillslope.def",
  zone = "defs/zone.def",
  soil = c("defs/soil_clay.def","defs/soil_clayloam.def","defs/soil_loam.def","defs/soil_loamysand.def","defs/soil_rock.def","defs/soil_sand.def","defs/soil_sandyclay.def","defs/soil_sandyclayloam.def","defs/soil_sandyloam.def","defs/soil_silt.def","defs/soil_siltyclay.def","defs/soil_siltyclayloam.def","defs/soil_siltyloam.def","defs/soil_water.def", "defs/soil_shallowloam.def", "defs/soil_shallowsandyclayloam.def", "defs/soil_shallowsandyloam.def"),
  landuse = c("defs/lu_undev.def","defs/lu_agriculture.def","defs/lu_raingarden.def","defs/lu_urban.def"),
  stratum = c("defs/veg_deciduous.def","defs/veg_deciduous_BES.def","defs/veg_eucalypt.def","defs/veg_evergreen.def","defs/veg_grass.def","defs/veg_lawn_2cm.def","defs/veg_lawn_5cm.def","defs/veg_lawn_10cm.def","defs/veg_nonveg.def"),
  basestations = "clim/cwtws32daymet.base")

Run RHESSys n times sequentially

Run RHESSys in Parallel using the foreach package, leaving 4 cores available for other tasks

Takes about 1 hour for 100 runs at 10x10m Daily WS32.

Read in RHESSys Calibration Outputs and Plot

Calibration is performed using GLUE methodology and Nash-Sutcliffe Model Efficiency (NSE), Log-Normal Nash-Sutcliffe Model Efficiency (lnNSE), and Kling-Gupta Efficiency (KGE).

setwd(system.file("extdata/", package = "RHESSysIOinR"))


# Read in RHESSys Calibration runs
for(i in 1:n) { assign(paste0("cwws32_run",i), readin_rhessys_output(paste0("out/cwws32_run",i)))}

#Prepare observed Streamflow data
obsws32$Date<- as.Date(obsws32$Date)

## plot observed data for period of record and calibration/validation periods as defined by SM data
{plot(obsws32$Date,obsws32$discharge_mm, type = 'l', xlim=c(as.Date(input_rhessys$start_date, "%Y%m%d"),as.Date(input_rhessys$end_date, "%Y%m%d")),ylab = "Observed Discharge mm/d", xlab = "Date", main = "Observed Streamflow Discharge for period of time model is run")
abline(v=as.Date(Caldates[1]), col = 'red', lwd = 2, lty = 2)
abline(v=as.Date(Caldates[2]), col = 'red', lwd = 2, lty = 2)
abline(v=as.Date(Valdates[1]), col = 'blue', lwd = 2, lty = 5)
abline(v=as.Date(Valdates[2]), col = 'blue', lwd = 2, lty = 5)
legend("topright",inset = 0.05, legend=c("Calibration","Validation"), col=c("red","blue"), lty=1)}

##merge calibration runs with observed so dates match up, instead of subsetting 
for(i in 1:n) { assign(paste0("simmerge",i), merge(obsws32,eval(parse(text = paste0("cwws32_run",i,"$bd"))), by.x = "Date", by.y="date", all = FALSE))
  assign(paste0("calsubset",i),eval(parse(text = paste0("simmerge",i)))[eval(parse(text = paste0("simmerge",i,"$Date"))) >= Caldates[1] & eval(parse(text = paste0("simmerge",i,"$Date"))) <= Caldates[2], ])}

for(i in 1:n) { assign(paste0("simmergesm",i), merge(obsws32sm,eval(parse(text = paste0("cwws32_run",i,"$bd"))), by.x = "Date", by.y="date", all = FALSE))
  assign(paste0("calsubsetsm",i),eval(parse(text = paste0("simmergesm",i)))[eval(parse(text = paste0("simmergesm",i,"$Date"))) >= Caldates[1] & eval(parse(text = paste0("simmergesm",i,"$Date"))) <= Caldates[2], ])}

# simmerge1<- merge(obsws32,eval(parse(text = paste0("cwws32_run",i,"$bd"))), by.x = "Date", by.y="date", all = FALSE)
#simmergesm1<- merge(observedsubsetsm,eval(parse(text = paste0("cwws32_run",i,"$bd"))), by.x = "Date", by.y="date", all = FALSE)


#subset observed data to calibration dates
# calsubset <- eval(parse(text = paste0("simmerge",i)))[eval(parse(text = paste0("simmerge",i,"$Date"))) >= Caldates[1] & eval(parse(text = paste0("simmerge",i,"$Date"))) <= Caldates[2], ]

#calsubsetsm <- observedsubsetsm[observedsubsetsm$Date >= Caldates[1] & observedsubsetsm$Date <= Caldates[2], ]
setwd(system.file("extdata/", package = "RHESSysIOinR"))


NSElist<- c()
lnNSElist<- c()
KGElist<- c()

smNSElist<- c()
smlnNSElist<- c()
smKGElist<- c()

for(i in 1:n) { assign(paste0("NSEobs",i), NSE(sim = as.numeric(eval(parse(text = paste0("simmerge",i,"$streamflow")))), obs = as.numeric(eval(parse(text = paste0("simmerge",i,"$observedstreamflow")))) ))
  NSElist[[i]]<- eval(parse(text = paste0("NSEobs",i))) }

for(i in 1:n) {  assign(paste0("lnNSEobs",i), NSE( sim = as.numeric(eval(parse(text = paste0("simmerge",i,"$streamflow")))), obs = as.numeric(eval(parse(text = paste0("simmerge",i,"$observedstreamflow")))), FUN = log, epsilon = "Pushpalatha2012", na.rm=TRUE))
  lnNSElist[[i]]<- eval(parse(text = paste0("lnNSEobs",i))) }

for(i in 1:n) {  assign(paste0("KGEobs",i), KGE(sim = as.numeric(eval(parse(text = paste0("simmerge",i,"$streamflow")))), obs = as.numeric(eval(parse(text = paste0("simmerge",i,"$observedstreamflow"))))))
  KGElist[[i]]<- eval(parse(text = paste0("KGEobs",i))) }

for(i in 1:n) {  assign(paste0("smNSEobs",i), NSE(sim = as.numeric(eval(parse(text = paste0("simmergesm",i,"$rz_storage","/simmergesm",i,"$rootdepth")))), obs = as.numeric(eval(parse(text = paste0("simmergesm",i,"$mergedsoilmoisture"))))))
  smNSElist[[i]]<- eval(parse(text = paste0("smNSEobs",i))) }

for(i in 1:n) {  assign(paste0("smlnNSEobs",i), NSE(sim = as.numeric(eval(parse(text = paste0("simmergesm",i,"$rz_storage","/simmergesm",i,"$rootdepth")))), obs = as.numeric(eval(parse(text = paste0("simmergesm",i,"$mergedsoilmoisture")))), FUN = log, epsilon = "Pushpalatha2012", na.rm=TRUE))
  smlnNSElist[[i]]<- eval(parse(text = paste0("smlnNSEobs",i))) }

for(i in 1:n) {  assign(paste0("smKGEobs",i), KGE(sim = as.numeric(eval(parse(text = paste0("simmergesm",i,"$rz_storage","/simmergesm",i,"$rootdepth")))), obs = as.numeric(eval(parse(text = paste0("simmergesm",i,"$mergedsoilmoisture"))))))
  smKGElist[[i]]<- eval(parse(text = paste0("smKGEobs",i))) }

NSElist[1:5]
## [[1]]
## [1] -0.003875825
## 
## [[2]]
## [1] 0.0767704
## 
## [[3]]
## [1] 0.1992447
## 
## [[4]]
## [1] -0.01049478
## 
## [[5]]
## [1] 0.4550887
KGElist[1:5]
## [[1]]
## [1] 0.4662955
## 
## [[2]]
## [1] 0.495591
## 
## [[3]]
## [1] 0.5451089
## 
## [[4]]
## [1] 0.4611981
## 
## [[5]]
## [1] 0.5980345
lnNSElist[1:5]
## [[1]]
## [1] 0.2647235
## 
## [[2]]
## [1] 0.2197867
## 
## [[3]]
## [1] 0.2514198
## 
## [[4]]
## [1] 0.1662837
## 
## [[5]]
## [1] 0.3505744
Reduce(max, NSElist)
## [1] 0.4876977
Reduce(max, lnNSElist)
## [1] 0.3570073
Reduce(max, KGElist)
## [1] 0.6270194
Reduce(max, smNSElist)
## [1] -1.624272
Reduce(max, smlnNSElist)
## [1] -1.50897
Reduce(max, smKGElist)
## [1] 0.002953977
which(NSElist == Reduce(max, NSElist), arr.ind=TRUE)
## [1] 142
which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE)
## [1] 58
which(KGElist == Reduce(max, KGElist), arr.ind=TRUE)
## [1] 65
which(smNSElist == Reduce(max, smNSElist), arr.ind=TRUE)
## [1] 72
which(smlnNSElist == Reduce(max, smlnNSElist), arr.ind=TRUE)
## [1] 34
which(smKGElist == Reduce(max, smKGElist), arr.ind=TRUE)
## [1] 4
paramtable<- read.table(file = "out/cwparamtable032823.txt")


paramtable[which(NSElist == Reduce(max, NSElist), arr.ind=TRUE),]
##             m        k  soil_dep      m_v      k_v      gw1       gw2
## 4947 1.478595 6.810474 0.3167354 1.225985 14.37205 0.266005 0.1388771
paramtable[which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE),]
##             m       k  soil_dep      m_v      k_v        gw1       gw2
## 6586 1.996897 11.9516 0.1659503 14.59814 11.96051 0.04187779 0.3730755
paramtable[which(KGElist == Reduce(max, KGElist), arr.ind=TRUE),]
##             m        k  soil_dep       m_v      k_v        gw1       gw2
## 5563 1.669432 26.16634 0.3261013 0.3147534 11.40159 0.03675462 0.2058475
paramtable[which(smNSElist == Reduce(max, smNSElist), arr.ind=TRUE),]
##             m        k  soil_dep       m_v      k_v       gw1       gw2
## 960 0.8611576 1.628167 0.4270699 0.5689802 12.58125 0.2944458 0.2843299
paramtable[which(smlnNSElist == Reduce(max, smlnNSElist), arr.ind=TRUE),]
##             m        k  soil_dep      m_v      k_v       gw1       gw2
## 6887 1.241347 2.344012 0.4243403 13.31544 6.071543 0.1391299 0.3108362
paramtable[which(smKGElist == Reduce(max, smKGElist), arr.ind=TRUE),]
##              m        k  soil_dep       m_v      k_v      gw1       gw2
## 9148 0.6367556 6.600494 0.2709346 0.7792454 11.35695 0.262354 0.4609033
# paramtable[7,]

###plot all data for GLUE

paramtable$NSE<- NSElist
paramtable$lnNSE<- lnNSElist
paramtable$KGE<- KGElist
paramtable$run<- seq.int(nrow(paramtable))

paramtable$smNSE<- smNSElist
paramtable$smlnNSE<- smlnNSElist
paramtable$smKGE<- smKGElist

paramtable<- as.data.frame(lapply(paramtable,unlist))


top100<- paramtable[order(-paramtable$lnNSE),][1:100,]

top100$run
##   [1]  58  19 126 125  57 137  40 117   5  76  16  94  42  50  79 128  80 127
##  [19] 146  51  49   8  62 142 108  92 135  26  66  41  78  69  87 101  91 111
##  [37]  45  73  88  27  71 133 136 144  99  89  21  28 107  67  22  81 110 122
##  [55]  20   7  68  85 114  33 145  15  25   6 119 118  53 140  59  98 121 143
##  [73]  44 123 120 124  30 116 113 109 139 102  65  47  83 134  14 104  82  29
##  [91] 112 100  64 131  77  48  63  84  36  52
#use values from top 100 run to table of daily streamflow predictions, for each day predict 95% prediction interval and mean then plot
#subset streamflow value for each day - day 1 day 2 etc
#put all values into a table

top100[["run"]]
##   [1]  58  19 126 125  57 137  40 117   5  76  16  94  42  50  79 128  80 127
##  [19] 146  51  49   8  62 142 108  92 135  26  66  41  78  69  87 101  91 111
##  [37]  45  73  88  27  71 133 136 144  99  89  21  28 107  67  22  81 110 122
##  [55]  20   7  68  85 114  33 145  15  25   6 119 118  53 140  59  98 121 143
##  [73]  44 123 120 124  30 116 113 109 139 102  65  47  83 134  14 104  82  29
##  [91] 112 100  64 131  77  48  63  84  36  52
toppredictions<- data.frame(cwws32_run1$bd$date)

for (i in top100$run){
 # print(i)
  toppredictions[,ncol(toppredictions)+1] <- eval(parse(text=paste0("cwws32_run",i,"$bd$streamflow")))
  colnames(toppredictions)[ncol(toppredictions)] <- paste0("stream",i)
}

head(toppredictions[1:5,])
##   cwws32_run1.bd.date stream58 stream19 stream126 stream125 stream57 stream137
## 1          2015-11-01 7.301086 7.621391  7.234245  7.033243 7.219342  6.935841
## 2          2015-11-02 8.339411 8.481753  8.223325  8.076945 8.454141  8.075405
## 3          2015-11-03 7.622130 7.759192  7.172434  7.837778 8.350942  8.026212
## 4          2015-11-04 7.194301 7.304438  6.680168  7.257826 7.691580  7.497977
## 5          2015-11-05 7.160050 7.311759  6.628511  7.158218 7.412052  7.315987
##   stream40 stream117  stream5 stream76 stream16 stream94 stream42 stream50
## 1 6.782780  7.173283 6.999768 7.718024 7.229691 7.243884 6.920345 7.543200
## 2 7.845375  8.105602 7.945108 8.758111 8.319140 8.427707 8.165312 8.586470
## 3 7.756190  7.604787 7.850628 7.978340 7.729866 8.074425 8.206078 7.850840
## 4 7.245324  7.235594 7.204286 7.431195 7.370403 7.595946 7.660926 7.434651
## 5 7.036828  7.255528 7.144896 7.330426 7.420635 7.441549 7.419535 7.375747
##   stream79 stream128 stream80 stream127 stream146 stream51 stream49  stream8
## 1 6.842465  7.031354 6.967658  6.908412  7.002993 6.982513 7.046788 7.197182
## 2 7.613180  7.979749 7.936525  7.709666  7.807816 7.912854 8.251385 8.231783
## 3 7.243560  7.685994 7.711997  7.469563  7.302992 7.773524 8.112560 7.743575
## 4 6.811250  7.303197 7.208213  6.940977  7.009744 7.280964 7.649492 7.392883
## 5 6.737169  7.224022 7.158457  6.960004  7.127518 7.109254 7.494514 7.438454
##   stream62 stream142 stream108 stream92 stream135 stream26 stream66 stream41
## 1 6.890182  7.480707  7.050082 6.994545  8.149884 6.128937 7.076216 7.746484
## 2 7.745624  8.359867  7.697038 8.156216  9.069572 7.083323 8.226110 8.935218
## 3 7.625534  7.664917  7.360530 7.866016  7.740992 6.885493 7.897024 8.100632
## 4 7.266316  7.213749  7.001604 7.509209  7.264247 6.467991 7.464590 7.517883
## 5 7.186921  7.162419  7.042933 7.440129  7.260217 6.335203 7.384197 7.443158
##   stream78 stream69 stream87 stream101 stream91 stream111 stream45  stream73
## 1 6.794378 7.119836 7.011195  6.720076 7.547154  6.723873 6.934614  9.046141
## 2 7.660612 7.907789 7.858544  7.854581 8.689763  7.517447 7.604283 10.119782
## 3 7.610500 7.633965 7.695322  7.900355 8.029957  6.980677 7.311125  9.198160
## 4 7.036635 7.211203 7.264788  7.409750 7.537334  6.705402 7.023489  8.598048
## 5 7.061996 7.198216 7.109192  7.301078 7.457512  6.901077 7.142527  8.370346
##   stream88 stream27 stream71 stream133 stream136 stream144 stream99 stream89
## 1 6.696809 7.066230 7.179298  6.772535  9.074359  6.780157 6.669808 6.877207
## 2 7.753083 8.544112 8.333969  7.501724  9.963068  7.873735 7.499996 7.675935
## 3 7.602220 8.365502 8.008033  6.912154  7.925216  7.414847 7.176466 7.560846
## 4 7.336603 7.767062 7.515260  6.728857  7.254993  7.135348 6.913612 7.235897
## 5 7.307955 7.521184 7.458419  6.866491  7.204661  7.105699 6.960972 7.206261
##   stream21 stream28 stream107 stream67 stream22 stream81 stream110 stream122
## 1 6.905437 6.645647  7.209404 6.768421 6.878164 6.894000  6.484788  7.172463
## 2 8.396006 7.596875  8.336647 7.498912 7.800314 7.756850  7.345957  8.224653
## 3 8.226518 7.715250  7.948836 7.313388 7.270891 7.206889  6.939228  7.731809
## 4 7.674125 7.161561  7.541502 7.062653 6.527299 7.006527  6.764322  7.364885
## 5 7.481645 7.145523  7.480826 7.150338 6.145675 7.169440  6.895039  7.450900
##   stream20  stream7 stream68 stream85 stream114 stream33 stream145 stream15
## 1 6.890570 7.051987 7.123813 6.690226  7.013109 7.260308  6.551573 6.539392
## 2 7.904469 8.205911 8.342977 7.342435  7.974642 8.258143  7.489137 7.057671
## 3 7.505268 7.950429 8.033743 7.126135  7.487044 7.591996  7.239523 6.481352
## 4 7.318275 7.511801 7.616806 6.868980  7.244088 7.344860  7.016527 6.276343
## 5 7.425878 7.432365 7.493687 7.006537  7.407971 7.456604  7.122548 6.412394
##   stream25  stream6 stream119 stream118 stream53 stream140 stream59 stream98
## 1 6.454897 6.731981  7.288448  6.697253 6.753820  6.710789 7.173417 6.622145
## 2 7.339007 7.870349  8.402888  7.489018 7.329145  7.718049 8.507835 7.299454
## 3 6.869175 8.070802  7.945091  7.234954 7.021659  7.934160 8.002741 6.976581
## 4 6.729247 7.459924  7.570510  6.919776 6.733350  7.424889 7.578688 6.747043
## 5 6.899262 7.264626  7.512899  6.909700 6.882012  7.272107 7.531734 6.941437
##   stream121 stream143 stream44 stream123 stream120 stream124 stream30 stream116
## 1  7.066062  7.158570 7.370372  8.067720  6.041898  7.276809 7.248389  6.801399
## 2  8.075929  8.437932 8.183350  9.316698  6.908755  8.509898 8.522186  7.648686
## 3  7.549531  8.219739 7.664192  8.133987  6.222903  7.870076 8.105675  7.320784
## 4  7.329247  7.653680 7.728082  7.512038  6.124654  7.566975 7.698694  6.945106
## 5  7.473999  7.550356 8.447302  7.424370  6.416691  7.614377 7.597161  6.954804
##   stream113 stream109 stream139 stream102 stream65 stream47 stream83 stream134
## 1  6.636456  6.606008  6.656626  7.891642 7.290588 7.798762 6.601900  7.746128
## 2  7.359144  7.289312  7.313785  9.290798 7.519508 9.089471 7.074206  8.950239
## 3  7.188012  7.038974  7.106471  8.240449 6.054104 8.014270 6.763529  8.003373
## 4  6.864055  6.837338  6.925474  7.678467 5.469140 7.481455 6.589236  7.595872
## 5  6.867512  6.936015  7.039053  7.584114 5.396264 7.433285 6.791125  7.641679
##   stream14 stream104 stream82 stream29 stream112 stream100 stream64 stream131
## 1 8.355276  7.406120 7.747981 7.056639  6.655416  7.285621 6.516596  7.940796
## 2 9.638516  8.769455 8.316809 7.520984  7.387885  8.541844 7.289585  9.125171
## 3 8.153606  8.202823 7.994366 6.672738  7.050887  8.056511 7.409187  7.844993
## 4 7.606671  7.699445 7.583742 6.337906  6.710635  7.659458 6.949297  7.451631
## 5 7.583419  7.620370 7.426955 6.275115  6.741510  7.637652 6.936601  7.540327
##   stream77 stream48 stream63 stream84 stream36  stream52
## 1 6.739532 6.564259 7.329269 7.255041 8.089297  9.153537
## 2 7.421621 7.321632 8.625799 8.632695 9.498303 10.266372
## 3 7.128002 7.116893 8.083832 8.221432 8.347417  7.932561
## 4 6.949546 6.892591 7.736889 7.786542 7.678440  7.309088
## 5 7.012069 6.889876 7.744266 7.698680 7.535520  7.382322
toppredictions<- toppredictions[,-1]

toppredictions$upper<- 0
toppredictions$mean<- 0
toppredictions$lower<- 0

for (i in 1:nrow(toppredictions)){
 #  print(i)
  upper= CI(as.numeric(toppredictions[i,]),ci=0.95)[1]
  mean= CI(as.numeric(toppredictions[i,]),ci=0.95)[2]
  lower= CI(as.numeric(toppredictions[i,]),ci=0.95)[3]
  
  toppredictions[i,]$upper<- upper
  toppredictions[i,]$mean<- mean
  toppredictions[i,]$lower<- lower
}

toppredictions<-cbind(date = cwws32_run1$bd$date,toppredictions)

head(toppredictions[1:5,])
##         date stream58 stream19 stream126 stream125 stream57 stream137 stream40
## 1 2015-11-01 7.301086 7.621391  7.234245  7.033243 7.219342  6.935841 6.782780
## 2 2015-11-02 8.339411 8.481753  8.223325  8.076945 8.454141  8.075405 7.845375
## 3 2015-11-03 7.622130 7.759192  7.172434  7.837778 8.350942  8.026212 7.756190
## 4 2015-11-04 7.194301 7.304438  6.680168  7.257826 7.691580  7.497977 7.245324
## 5 2015-11-05 7.160050 7.311759  6.628511  7.158218 7.412052  7.315987 7.036828
##   stream117  stream5 stream76 stream16 stream94 stream42 stream50 stream79
## 1  7.173283 6.999768 7.718024 7.229691 7.243884 6.920345 7.543200 6.842465
## 2  8.105602 7.945108 8.758111 8.319140 8.427707 8.165312 8.586470 7.613180
## 3  7.604787 7.850628 7.978340 7.729866 8.074425 8.206078 7.850840 7.243560
## 4  7.235594 7.204286 7.431195 7.370403 7.595946 7.660926 7.434651 6.811250
## 5  7.255528 7.144896 7.330426 7.420635 7.441549 7.419535 7.375747 6.737169
##   stream128 stream80 stream127 stream146 stream51 stream49  stream8 stream62
## 1  7.031354 6.967658  6.908412  7.002993 6.982513 7.046788 7.197182 6.890182
## 2  7.979749 7.936525  7.709666  7.807816 7.912854 8.251385 8.231783 7.745624
## 3  7.685994 7.711997  7.469563  7.302992 7.773524 8.112560 7.743575 7.625534
## 4  7.303197 7.208213  6.940977  7.009744 7.280964 7.649492 7.392883 7.266316
## 5  7.224022 7.158457  6.960004  7.127518 7.109254 7.494514 7.438454 7.186921
##   stream142 stream108 stream92 stream135 stream26 stream66 stream41 stream78
## 1  7.480707  7.050082 6.994545  8.149884 6.128937 7.076216 7.746484 6.794378
## 2  8.359867  7.697038 8.156216  9.069572 7.083323 8.226110 8.935218 7.660612
## 3  7.664917  7.360530 7.866016  7.740992 6.885493 7.897024 8.100632 7.610500
## 4  7.213749  7.001604 7.509209  7.264247 6.467991 7.464590 7.517883 7.036635
## 5  7.162419  7.042933 7.440129  7.260217 6.335203 7.384197 7.443158 7.061996
##   stream69 stream87 stream101 stream91 stream111 stream45  stream73 stream88
## 1 7.119836 7.011195  6.720076 7.547154  6.723873 6.934614  9.046141 6.696809
## 2 7.907789 7.858544  7.854581 8.689763  7.517447 7.604283 10.119782 7.753083
## 3 7.633965 7.695322  7.900355 8.029957  6.980677 7.311125  9.198160 7.602220
## 4 7.211203 7.264788  7.409750 7.537334  6.705402 7.023489  8.598048 7.336603
## 5 7.198216 7.109192  7.301078 7.457512  6.901077 7.142527  8.370346 7.307955
##   stream27 stream71 stream133 stream136 stream144 stream99 stream89 stream21
## 1 7.066230 7.179298  6.772535  9.074359  6.780157 6.669808 6.877207 6.905437
## 2 8.544112 8.333969  7.501724  9.963068  7.873735 7.499996 7.675935 8.396006
## 3 8.365502 8.008033  6.912154  7.925216  7.414847 7.176466 7.560846 8.226518
## 4 7.767062 7.515260  6.728857  7.254993  7.135348 6.913612 7.235897 7.674125
## 5 7.521184 7.458419  6.866491  7.204661  7.105699 6.960972 7.206261 7.481645
##   stream28 stream107 stream67 stream22 stream81 stream110 stream122 stream20
## 1 6.645647  7.209404 6.768421 6.878164 6.894000  6.484788  7.172463 6.890570
## 2 7.596875  8.336647 7.498912 7.800314 7.756850  7.345957  8.224653 7.904469
## 3 7.715250  7.948836 7.313388 7.270891 7.206889  6.939228  7.731809 7.505268
## 4 7.161561  7.541502 7.062653 6.527299 7.006527  6.764322  7.364885 7.318275
## 5 7.145523  7.480826 7.150338 6.145675 7.169440  6.895039  7.450900 7.425878
##    stream7 stream68 stream85 stream114 stream33 stream145 stream15 stream25
## 1 7.051987 7.123813 6.690226  7.013109 7.260308  6.551573 6.539392 6.454897
## 2 8.205911 8.342977 7.342435  7.974642 8.258143  7.489137 7.057671 7.339007
## 3 7.950429 8.033743 7.126135  7.487044 7.591996  7.239523 6.481352 6.869175
## 4 7.511801 7.616806 6.868980  7.244088 7.344860  7.016527 6.276343 6.729247
## 5 7.432365 7.493687 7.006537  7.407971 7.456604  7.122548 6.412394 6.899262
##    stream6 stream119 stream118 stream53 stream140 stream59 stream98 stream121
## 1 6.731981  7.288448  6.697253 6.753820  6.710789 7.173417 6.622145  7.066062
## 2 7.870349  8.402888  7.489018 7.329145  7.718049 8.507835 7.299454  8.075929
## 3 8.070802  7.945091  7.234954 7.021659  7.934160 8.002741 6.976581  7.549531
## 4 7.459924  7.570510  6.919776 6.733350  7.424889 7.578688 6.747043  7.329247
## 5 7.264626  7.512899  6.909700 6.882012  7.272107 7.531734 6.941437  7.473999
##   stream143 stream44 stream123 stream120 stream124 stream30 stream116 stream113
## 1  7.158570 7.370372  8.067720  6.041898  7.276809 7.248389  6.801399  6.636456
## 2  8.437932 8.183350  9.316698  6.908755  8.509898 8.522186  7.648686  7.359144
## 3  8.219739 7.664192  8.133987  6.222903  7.870076 8.105675  7.320784  7.188012
## 4  7.653680 7.728082  7.512038  6.124654  7.566975 7.698694  6.945106  6.864055
## 5  7.550356 8.447302  7.424370  6.416691  7.614377 7.597161  6.954804  6.867512
##   stream109 stream139 stream102 stream65 stream47 stream83 stream134 stream14
## 1  6.606008  6.656626  7.891642 7.290588 7.798762 6.601900  7.746128 8.355276
## 2  7.289312  7.313785  9.290798 7.519508 9.089471 7.074206  8.950239 9.638516
## 3  7.038974  7.106471  8.240449 6.054104 8.014270 6.763529  8.003373 8.153606
## 4  6.837338  6.925474  7.678467 5.469140 7.481455 6.589236  7.595872 7.606671
## 5  6.936015  7.039053  7.584114 5.396264 7.433285 6.791125  7.641679 7.583419
##   stream104 stream82 stream29 stream112 stream100 stream64 stream131 stream77
## 1  7.406120 7.747981 7.056639  6.655416  7.285621 6.516596  7.940796 6.739532
## 2  8.769455 8.316809 7.520984  7.387885  8.541844 7.289585  9.125171 7.421621
## 3  8.202823 7.994366 6.672738  7.050887  8.056511 7.409187  7.844993 7.128002
## 4  7.699445 7.583742 6.337906  6.710635  7.659458 6.949297  7.451631 6.949546
## 5  7.620370 7.426955 6.275115  6.741510  7.637652 6.936601  7.540327 7.012069
##   stream48 stream63 stream84 stream36  stream52    upper     mean    lower
## 1 6.564259 7.329269 7.255041 8.089297  9.153537 7.173589 6.915854 6.658119
## 2 7.321632 8.625799 8.632695 9.498303 10.266372 8.163679 7.866137 7.568594
## 3 7.116893 8.083832 8.221432 8.347417  7.932561 7.676408 7.406455 7.136502
## 4 6.892591 7.736889 7.786542 7.678440  7.309088 7.268135 7.015501 6.762866
## 5 6.889876 7.744266 7.698680 7.535520  7.382322 7.252374 7.001687 6.751001
# toppredictions[1:500,]

ggplot(toppredictions,aes(date,mean,group=1))+
  geom_line(col = 'red')+
  geom_ribbon(aes(ymin=lower, ymax=upper), fill = "blue", alpha=0.3)+
  geom_line(data=calsubset1,aes(x=Date,y=discharge_mm))+
  ylab("Discharge mm")+
  xlab("Date")+
  labs(title = "Top 100 Model Runs vs Measured Streamflow")

par(mfrow=c(1, 1))

{plot(paramtable$m,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - M", main = "M - Decay of hydraulic conductivity with depth")
abline(h=0,lty=2)}

{plot(paramtable$k,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - K", main = "K - Hydraulic conductivity at the surface")
abline(h=0,lty=2)}

{plot(paramtable$soil_dep,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - Sd", main = "Sd - soil depth")
abline(h=0,lty=2)}

{plot(paramtable$m_v,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - Mv", main = "Mv - Vertical decay of hydraulic conductivity with depth")
abline(h=0,lty=2)}

{plot(paramtable$k_v,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - Kv", main = "Kv - Vertical hydraulic conductivity at the surface")
abline(h=0,lty=2)}

{plot(paramtable$gw1,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - GW1", main = "GW1 - Saturated to Groundwater Coefficient")
abline(h=0,lty=2)}

{plot(paramtable$gw2,paramtable$NSE, pch = 20,ylab="Nash-Sutcliffe model efficiency coefficient", xlab = "Multiplier - GW2", main = "GW2 - Groundwater Loss Coefficient")
abline(h=0,lty=2)}

old.par <- par(mfrow=c(3, 7), mar = c(5,4,4,2))

plot(paramtable$m,paramtable$NSE,xlab = "m parameter", pch = 20)
lines(lowess(paramtable$m,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))
plot(paramtable$k,paramtable$NSE,xlab = "lateral ksat parameter", pch = 20)
lines(lowess(paramtable$k,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))
plot(paramtable$soil_dep,paramtable$NSE,xlab = "sdepth parameter", pch = 20)
lines(lowess(paramtable$soil_dep,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))
plot(paramtable$m_v,paramtable$NSE,xlab = "mv parameter", pch = 20)
lines(lowess(paramtable$m_v,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))
plot(paramtable$k_v,paramtable$NSE,xlab = "ksatv parameter", pch = 20)
lines(lowess(paramtable$k_v,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))
plot(paramtable$gw1,paramtable$NSE,xlab = "gw1 parameter", pch = 20)
lines(lowess(paramtable$gw1,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))
plot(paramtable$gw2,paramtable$NSE,xlab = "gw2 parameter", pch = 20)
lines(lowess(paramtable$gw2,paramtable$NSE, f = 0.1),col="red",ylim = c(-1,0.7))


plot(paramtable$m,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "m parameter", pch = 20)
lines(lowess(paramtable$m,paramtable$lnNSE, f = 0.1),col="red")
plot(paramtable$k,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "lateral ksat parameter", pch = 20)
lines(lowess(paramtable$k,paramtable$lnNSE, f = 0.1),col="red")
plot(paramtable$soil_dep,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "sdepth parameter", pch = 20)
lines(lowess(paramtable$soil_dep,paramtable$lnNSE, f = 0.1),col="red")
plot(paramtable$m_v,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "mv parameter", pch = 20)
lines(lowess(paramtable$m_v,paramtable$lnNSE, f = 0.1),col="red")
plot(paramtable$k_v,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "ksatv parameter", pch = 20)
lines(lowess(paramtable$k_v,paramtable$lnNSE, f = 0.1),col="red")
plot(paramtable$gw1,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "gw1 parameter", pch = 20)
lines(lowess(paramtable$gw1,paramtable$lnNSE, f = 0.1),col="red")
plot(paramtable$gw2,paramtable$lnNSE,ylim = c(-1,0.7),xlab = "gw2 parameter", pch = 20)
lines(lowess(paramtable$gw2,paramtable$lnNSE, f = 0.1),col="red")


plot(paramtable$m,paramtable$KGE,ylim = c(-0.25,1),xlab = "m parameter", pch = 20)
lines(lowess(paramtable$m,paramtable$KGE, f = 0.1),col="red")
plot(paramtable$k,paramtable$KGE,ylim = c(-0.25,1),xlab = "lateral ksat parameter", pch = 20)
lines(lowess(paramtable$k,paramtable$KGE, f = 0.1),col="red")
plot(paramtable$soil_dep,paramtable$KGE,ylim = c(-0.25,1),xlab = "sdepth parameter", pch = 20)
lines(lowess(paramtable$soil_dep,paramtable$KGE, f = 0.1),col="red")
plot(paramtable$m_v,paramtable$KGE,ylim = c(-0.25,0.7),xlab = "mv parameter", pch = 20)
lines(lowess(paramtable$m_v,paramtable$KGE, f = 0.1),col="red")
plot(paramtable$k_v,paramtable$KGE,ylim = c(-0.25,0.7),xlab = "ksatv parameter", pch = 20)
lines(lowess(paramtable$k_v,paramtable$KGE, f = 0.1),col="red")
plot(paramtable$gw1,paramtable$KGE,ylim = c(-0.25,0.7),xlab = "gw1 parameter", pch = 20)
lines(lowess(paramtable$gw1,paramtable$KGE, f = 0.1),col="red")
plot(paramtable$gw2,paramtable$KGE,ylim = c(-0.25,0.7),xlab = "gw2 parameter", pch = 20)
lines(lowess(paramtable$gw2,paramtable$KGE, f = 0.1),col="red")

par(old.par)

old.par <- par(mfrow=c(1, 7), mar = c(5,4,4,2))

plot(paramtable$m,paramtable$smKGE,ylim = c(-0.25,1),xlab = "m parameter", pch = 20)
lines(lowess(paramtable$m,paramtable$smKGE, f = 0.1),col="red")
plot(paramtable$k,paramtable$smKGE,ylim = c(-0.25,1),xlab = "lateral ksat parameter", pch = 20)
lines(lowess(paramtable$k,paramtable$smKGE, f = 0.1),col="red")
plot(paramtable$soil_dep,paramtable$smKGE,ylim = c(-0.25,1),xlab = "sdepth parameter", pch = 20)
lines(lowess(paramtable$soil_dep,paramtable$smKGE, f = 0.1),col="red")
plot(paramtable$m_v,paramtable$smKGE,ylim = c(-0.25,0.7),xlab = "mv parameter", pch = 20)
lines(lowess(paramtable$m_v,paramtable$smKGE, f = 0.1),col="red")
plot(paramtable$k_v,paramtable$smKGE,ylim = c(-0.25,0.7),xlab = "ksatv parameter", pch = 20)
lines(lowess(paramtable$k_v,paramtable$smKGE, f = 0.1),col="red")
plot(paramtable$gw1,paramtable$smKGE,ylim = c(-0.25,0.7),xlab = "gw1 parameter", pch = 20)
lines(lowess(paramtable$gw1,paramtable$smKGE, f = 0.1),col="red")
plot(paramtable$gw2,paramtable$smKGE,ylim = c(-0.25,0.7),xlab = "gw2 parameter", pch = 20)
lines(lowess(paramtable$gw2,paramtable$smKGE, f = 0.1),col="red")

par(old.par)

### end plot all data for GLUE

paste0("cwws32_run",which(NSElist == Reduce(max, NSElist), arr.ind=TRUE),"$bd$streamflow")
## [1] "cwws32_run142$bd$streamflow"
paste0("cwws32_run",which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE),"$bd$streamflow")
## [1] "cwws32_run58$bd$streamflow"
paste0("cwws32_run",which(KGElist == Reduce(max, KGElist), arr.ind=TRUE),"$bd$streamflow")
## [1] "cwws32_run65$bd$streamflow"
{plot(cwws32_run33$bd$date,cwws32_run33$bd$streamflow, type = "l", col = 'red', main = "Observed Streamflow vs Top Runs", xlab = "Date", ylab = "Streamflow")
lines(cwws32_run65$bd$date,cwws32_run65$bd$streamflow, col = 'blue')
#lines(cwws32_run46$bd$date,cwws32_run46$bd$streamflow, col = 'green')
lines(simmerge1$Date,simmerge1$observedstreamflow, col = 'orange', lwd = 2, lty = 2)}

#{plot(paste("cwws32_run",which(NSElist == Reduce(max, NSElist), arr.ind=TRUE),"$bd$date", sep = ""),paste("cwws32_run",which(NSElist == Reduce(max, NSElist), arr.ind=TRUE),"$bd$streamflow", sep = ""), type = "l", col = 'red', main = "Observed Streamflow vs Top Runs", xlab = "Date", ylab = "Streamflow")
#lines(paste("cwws32_run",which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE),"$bd$date", sep = ""),paste("cwws32_run",which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE),"$bd$streamflow", sep = ""), col = 'blue')
#lines(cwws32_run80$bd$date,cwws32_run80$bd$streamflow, col = 'green')
#lines(paste("cwws32_run",which(KGElist == Reduce(max, KGElist), arr.ind=TRUE),"$bd$date", sep = ""),paste("cwws32_run",which(KGElist == Reduce(max, KGElist), arr.ind=TRUE),"$bd$streamflow", sep = ""), col = 'orange', lwd = 2, lty = 2)
#legend("topleft",inset = 0.05, legend=c("Oberved Discharge","Best NSE","BEST lnNSE","Best KGE"), col=c("orange","red","blue","green"), lty=c(2,1,1,1))}

{plot(cwws32_run13$bd$date,cwws32_run13$bd$streamflow, type = "l", col = 'red', main = "Observed Streamflow vs Top Runs", xlab = "Date", ylab = "Streamflow", ylim = c(0,10))
#lines(cwws32_run47$bd$date,cwws32_run47$bd$streamflow, col = 'blue')
lines(obsws32$Date,obsws32$discharge_mm, col = 'orange', lwd = 2, lty = 2)}

*** VALIDATION ***

Run Model once for validation, copied from previous code must clean up

Configure RHESSys Inputs/Outputs for a single run.

Use World Statefile created during spin-up. To run model one time at patch scale.

COMMAND LINE OPTIONS

-b Basin output option. Print out response variables for specified basins. -c Canopy stratum output option. Print out response variables for specified strata. -g Grow option. Try to read in dynamic bgc input data and output dynamic bgc parameters. -h Hillslope output option. Print out response variables for specified hillslopes. -p Patch output option. Print out response variables for specified patches. -r Routing option. Gives name of flow_table to define explicit routing connectivity. Also trigger use of explicit routing over TOPMODEL approach. -c Stratum output option. Print out response variables for specified strata.

input_rhessys = IOin_rhessys_input(
  version = rh_path,
  tec_file = "tecfiles/tec_daily",
  world_file = "CWWS32.world.Y2018M10D31H1.state",
  world_hdr_prefix = "CWWS32",
  flowtable = "CWWS321.flow",
  start = "2014 11 1 1",
  end = "2018 11 1 1",
  output_folder = "out",
  output_prefix = "cwws32",
  commandline_options = c("-b -g -p"))
 
# commandline_options = c("-b -g -p 1 128 154914 154914")

## TEC file dictates model output, begin output a year in to allow model SM to stabilize
# do not output_state or worldfile may be overwritten as output is created
input_tec_data = IOin_tec_std(start = "2015 11 1 1",
                              end = "2018 11 1 1",
                              output_state = FALSE)

input_hdr = IOin_hdr(
  basin = "defs/basin.def",
  hillslope = "defs/hillslope.def",
  zone = "defs/zone.def",
  soil = c("defs/soil_clay.def","defs/soil_clayloam.def","defs/soil_loam.def","defs/soil_loamysand.def","defs/soil_rock.def","defs/soil_sand.def","defs/soil_sandyclay.def","defs/soil_sandyclayloam.def","defs/soil_sandyloam.def","defs/soil_silt.def","defs/soil_siltyclay.def","defs/soil_siltyclayloam.def","defs/soil_siltyloam.def","defs/soil_water.def", "defs/soil_shallowloam.def", "defs/soil_shallowsandyclayloam.def", "defs/soil_shallowsandyloam.def"),
  landuse = "defs/lu_undev.def",
  stratum = "defs/veg_deciduous.def",
  basestations = "clim/cwtws32daymet.base")

## Calibrated Values 


which(NSElist == Reduce(max, NSElist), arr.ind=TRUE)
## [1] 142
which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE)
## [1] 58
which(KGElist == Reduce(max, KGElist), arr.ind=TRUE)
## [1] 65
which(smNSElist == Reduce(max, smNSElist), arr.ind=TRUE)
## [1] 72
which(smlnNSElist == Reduce(max, smlnNSElist), arr.ind=TRUE)
## [1] 34
which(smKGElist == Reduce(max, smKGElist), arr.ind=TRUE)
## [1] 4
paramtable[which(NSElist == Reduce(max, NSElist), arr.ind=TRUE),]
##            m        k  soil_dep      m_v      k_v      gw1       gw2       NSE
## 142 1.478595 6.810474 0.3167354 1.225985 14.37205 0.266005 0.1388771 0.4876977
##         lnNSE       KGE run     smNSE   smlnNSE      smKGE
## 142 0.3447955 0.5987777 142 -4.221292 -2.951442 -0.2579614
paramtable[which(lnNSElist == Reduce(max, lnNSElist), arr.ind=TRUE),]
##           m       k  soil_dep      m_v      k_v        gw1       gw2       NSE
## 58 1.996897 11.9516 0.1659503 14.59814 11.96051 0.04187779 0.3730755 0.4629913
##        lnNSE       KGE run     smNSE   smlnNSE      smKGE
## 58 0.3570073 0.6166828  58 -5.206479 -3.536803 -0.4590018
paramtable[which(KGElist == Reduce(max, KGElist), arr.ind=TRUE),]
##           m        k  soil_dep       m_v      k_v        gw1       gw2
## 65 1.669432 26.16634 0.3261013 0.3147534 11.40159 0.03675462 0.2058475
##          NSE     lnNSE       KGE run     smNSE   smlnNSE      smKGE
## 65 0.4219818 0.3230267 0.6270194  65 -4.489804 -3.079227 -0.1218914
paramtable[which(smNSElist == Reduce(max, smNSElist), arr.ind=TRUE),]
##            m        k  soil_dep       m_v      k_v       gw1       gw2
## 72 0.8611576 1.628167 0.4270699 0.5689802 12.58125 0.2944458 0.2843299
##            NSE     lnNSE       KGE run     smNSE   smlnNSE       smKGE
## 72 -0.04066875 0.2492639 0.4594109  72 -1.624272 -2.020654 -0.05025428
paramtable[which(smlnNSElist == Reduce(max, smlnNSElist), arr.ind=TRUE),]
##           m        k  soil_dep      m_v      k_v       gw1       gw2        NSE
## 34 1.241347 2.344012 0.4243403 13.31544 6.071543 0.1391299 0.3108362 0.09628193
##        lnNSE       KGE run     smNSE  smlnNSE      smKGE
## 34 0.2993379 0.5115521  34 -1.887176 -1.50897 -0.1419453
paramtable[which(smKGElist == Reduce(max, smKGElist), arr.ind=TRUE),]
##           m        k  soil_dep       m_v      k_v      gw1       gw2
## 4 0.6367556 6.600494 0.2709346 0.7792454 11.35695 0.262354 0.4609033
##           NSE     lnNSE       KGE run    smNSE   smlnNSE       smKGE
## 4 -0.01049478 0.1662837 0.4611981   4 -1.69059 -1.517735 0.002953977
stdpars<- IOin_std_pars(
            m = 0.64,
            k = 6.6,
            soil_dep= 0.27,
            m_v = 0.78,
            k_v = 11.36,
            gw1 = 0.26,
            gw2 = 0.46)

Run RHESSys once for validation

Read RHESSys Single Run Output

getwd()
## [1] "C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata"
singlerunvalidation<- readin_rhessys_output("out/cwws32_runvalidation")

plot(singlerunvalidation$bd$streamflow~singlerunvalidation$bd$date, type = "l", main = "Streamflow", xlab = "Date", ylab = "Streamflow", col = 'DarkBlue')

plot(as.Date(singlerunvalidation$bd$date),singlerunvalidation$bd$rz_storage, type = "l", col = 'black', main = "Basin Scale Root Zone Storage",
    xlab = "Date", ylab = "mm")

#basin scale soil moisture
plot(singlerunvalidation$bd$rz_storage/singlerunvalidation$bdg$root_depth~singlerunvalidation$bd$date, type = "l", main = "RZ_Storage/Root_Depth x Date", xlab = "Date", ylab = "rz_Storage/root_depth", col = 'brown')

plot(singlerunvalidation$bd$lai~singlerunvalidation$bd$date, type = "l", main = "LAI", xlab = "Date", ylab = "LAI", col = 'DarkGreen')

plot(singlerunvalidation$bd$pet~singlerunvalidation$bd$date, type = "l", main = "Potential Evapotranspiration", xlab = "Date", ylab = "PET", col = 'DarkBlue')

plot(singlerunvalidation$bd$et~singlerunvalidation$bd$date, type = "l", main = "Evapotranspiration", xlab = "Date", ylab = "ET", col = 'darkslategray')

plot((singlerunvalidation$bd$unsat_stor/singlerunvalidation$bd$sat_def_z)~singlerunvalidation$bd$date, type = "l", main = "unsat_stor/sat_def_z", xlab = "Date", ylab = "vwc", col = 'DarkGreen')

Merge Observed and Simulated Data

## merge values instead of subsetting to match up with simulated dates
singlerunvalidation$bd<- merge(singlerunvalidation$bd,obsws32, by.x = "date", by.y = "Date", all = FALSE)

## merge soil moisture values but include all data
singlerunvalidation$data <-merge(singlerunvalidation$bd,obsws32sm, by.x = "date", by.y = "Date", all = TRUE)

Plot Observed and Simulated Data, Plot RHESSys Outputs

## Plot Hydrograph comparing modeled and observed streamflow
{hydrograph(input=singlerunvalidation$bd, streamflow=singlerunvalidation$bd$observedstreamflow, streamflow2 = singlerunvalidation$bd$streamflow, timeSeries = singlerunvalidation$bd$date, precip = singlerunvalidation$bd$streamflow,
           P.units = "mm", S.units = "mm normalized by basin area", S1.col = 'Blue', S2.col = 'Red')

legend("topleft",inset = .2, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("blue","red"), lty=1:2)}

{plot(singlerunvalidation$data$date,singlerunvalidation$data$observedstreamflow.x, type = "l", xlab = "Date", ylab = "Streamflow", main = "Predicted vs Observed Streamflow")
lines(singlerunvalidation$data$date, singlerunvalidation$data$streamflow, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:1)}

{plot(singlerunvalidation$data$date,singlerunvalidation$data$observedstreamflow.x, type = "l", xlab = "Date", ylab = "Streamflow", main = "Predicted vs Observed Log Streamflow", log = "y")
lines(singlerunvalidation$data$date, singlerunvalidation$data$streamflow, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:1)}

predictedvsobservedstreamlm<- lm(observedstreamflow~streamflow, data = singlerunvalidation$bd)

{plot(singlerunvalidation$bd$observedstreamflow,singlerunvalidation$bd$streamflow, xlab = "Observed", ylab = "Predicted", main = "Predicted vs Observed Streamflow", xlim = c(0,55), ylim = c(0,55))
abline(a=0, b=1, col = 'RED', lty = 2, lwd = 2)
legend("right",inset = 0.01, legend = paste("R2 =",format(summary(predictedvsobservedstreamlm)$r.squared,digits=3)))}

{plot(singlerunvalidation$data$date,(singlerunvalidation$data$rz_storage/singlerunvalidation$data$rootdepth), type = "l", ylim = c(0,0.35), xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Soil Moisture")
lines(singlerunvalidation$data$date, singlerunvalidation$data$mergedsoilmoisture, col = 'RED')
legend("topright",inset = 0, legend=c("Observed Soil Moisture","Predicted Soil Moisture"), col=c("black","red"), lty=1:1)}

predictedvsobservedsoillm<- lm((rz_storage/rootdepth)~mergedsoilmoisture, data = singlerunvalidation$data)

{plot((singlerunvalidation$data$rz_storage/singlerunvalidation$data$rootdepth),singlerunvalidation$data$mergedsoilmoisture, xlab = "Observed", ylab = "Predicted", main = "Predicted vs Observed Basin Scale Soil Moisture", xlim = c(0.05,0.4), ylim = c(0.05,0.4))
abline(a=0, b=1, col = 'RED', lty = 2, lwd = 2)
legend("right",inset = 0.01, legend = paste("R2 =",format(summary(predictedvsobservedsoillm)$r.squared,digits=3)))}

## Plot other model outputs

par(mfrow=c(5,1), mar = c(2,5,4,2))

plot(singlerunvalidation$bd$date,singlerunvalidation$bd$tmin, type = "l", col = 'BLUE', ylim = c(-15,30), ylab = "Temperature", xlab = "Date")
lines(singlerunvalidation$bd$date,singlerunvalidation$bd$tmax, type = "l", col = 'RED')
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$precip, type = "h", ylab = "Precipitation", xlab = "Date")
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$streamflow, type = "l", ylab = "Streamflow", xlab = "Date", col = 'blue')
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$baseflow, type = "l", ylab = "Baseflow", xlab = "Date", col = 'brown')
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$rz_storage, type = "l", ylab = "RZ Storage", xlab = "Date", col = 'dark green')

plot(singlerunvalidation$bd$date,singlerunvalidation$bd$psn ,type = "l", ylab = "PSN", xlab = "Date")
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$et, type = "l", ylab = "Evaoptranspiration", xlab = "Date")
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$pet, type = "l", ylab = "PET", xlab = "Date")
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$gw.Qout, type = "l", ylab = "Q out", xlab = "Date")
plot(singlerunvalidation$bd$date,singlerunvalidation$bd$gw.storage, type = "l", ylab = "GW Storage", xlab = "Date")

par(mfrow=c(1,1))

### will only plot in growth mode , used to check that vegetation is initialized properly ##
plot(singlerunvalidation$bdg$date, singlerunvalidation$bdg$lai, type = "l", main = "LAI", xlab = "Date", col = 'dark green')

plot(singlerunvalidation$bdg$date, singlerunvalidation$bdg$soilc, type = "l", main = "Soil Carbon", xlab = "Date", col = 'brown')

## plot cal/val sets

{plot(singlerunvalidation$data$date,(singlerunvalidation$data$rz_storage/singlerunvalidation$data$rootdepth), type = "l", ylim = c(0,0.35), xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Soil Moisture")
lines(obsws32smcal$Date,obsws32smcal$mergedsoilmoisture, col = "BLUE")
lines(obsws32smval$Date,obsws32smval$mergedsoilmoisture, col = "RED")
legend("topright",inset = 0, legend=c("Predicted Soil Moisture","Observed Calibration Soil Moisture", "Observed Validation Soil Moisture"), col=c("black","blue","red"), lty=c(1,1,1))}

{plot(singlerunvalidation$data$date,singlerunvalidation$data$streamflow, type = "l", xlab = "Date", ylab = "RZ_Storage/Rootdepth", main = "Predicted vs Observed Basin Scale Streamflow")
lines(obsws32cal$Date,obsws32cal$observedstreamflow, col = "BLUE")
lines(obsws32val$Date,obsws32val$observedstreamflow, col = "RED")
legend("topright",inset = 0, legend=c("Predicted Streamflow","Observed Calibration Streamflow", "Observed Validation Streamflow"), col=c("black","blue","red"), lty=c(1,1,1))}

Spatial Data import and Display

# Original input maps have been cropped before processing in R, without cropping the RHESSYs preparation in the R environment will not function properly, below are the input maps used by R
ws32<- raster("grassexport/cwt_ws32/basin_ws32.tif")
patches<- raster("grassexport/cwt_ws32/patch.tif")
acc10<- raster("grassexport/cwt_ws32/acc10.tif")
aspect10<- raster("grassexport/cwt_ws32/aspect10.tif")
#bt1000<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/b.t1000.tif")
dem10<- raster("grassexport/cwt_ws32/dem10.tif")
dem10f<- raster("grassexport/cwt_ws32/dem10f.tif")
dir10<- raster("grassexport/cwt_ws32/dir10.tif")
drain10<- raster("grassexport/cwt_ws32/drain10.tif")
#ht1000<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/h.t1000.tif")
hillslope<- raster("grassexport/cwt_ws32/hillslope.tif")
roads<- brick("grassexport/cwt_ws32/roads.tif")
#shade10<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/shade10.tif")
slope10<- raster("grassexport/cwt_ws32/slope10.tif")
streams<- raster("grassexport/cwt_ws32/streams.tif")
topidx10<- raster("grassexport/cwt_ws32/topidx10_100.tif")
#xmap<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/xmap.tif")
#ymap<- raster("C:/Users/Carlos/Documents/R/win-library/4.0/RHESSysIOinR/extdata/rasters/cwt_ws32/ws32/ymap.tif")


#plot(ws32, main ="ws32")
#plot(patches, main = "patches")
#plot(acc10, main = "acc10")
#plot(aspect10, main = "aspect10")
#levelplot(bt1000, col.regions = brewer.pal(name='Spectral', n = 11),at=seq(0,6500,1))
#plot(bt1000, main = "bt1000")
#plot(dem10, main = "dem10")
#plot(dem10f, main = "dem10f")
#plot(dir10, main = "dir10")
#plot(drain10, main = "drain10")
#plot(ht1000, main = "ht1000")
#plot(hillslope, main = "hillslope")
#plot(roads, main = "roads")
#plot(shade10, main = "shade10")
#plot(slope10, main = "slope10")
#plot(streams, main = "streams")
#plot(topidx10, main = "topidx10")
#plot(xmap, main = "xmap")
#plot(ymap, main = "ymap")
{plot(dem10, col = colorramp, main = "Coweeta Observed Soil Moisture Locations")
points(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1]), c(ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", lwd = 1)}

matrix(c(ObservedSMLocations$Name[3],ObservedSMLocations$Name[4],ObservedSMLocations$Name[5],ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), nrow = 3, ncol = 3)
##      [,1]              [,2]               [,3]              
## [1,] "132 plot center" "275432.649744231" "3881522.05364127"
## [2,] "232 plot center" "275578.214010129" "3881428.55011708"
## [3,] "332 plot center" "275760.864807773" "3881284.18476691"
## Extract values for single date
displayday <- singlerunvalidation$pd[which(singlerunvalidation$pd$date=="2018-03-21")]

## reclassify all patches to na and then reclassify to model output values
reclass_all<- c(-2147483648, 2147483647,NA)
reclass_allm<- matrix(reclass_all, ncol = 3, byrow = TRUE)
displaydaypatches<- reclassify(patches,reclass_allm)

displayday$calculatedrzsm <- displayday$rz_storage/displayday$root.depth

reclass_displayday <- cbind(displayday$patchID,displayday$calculatedrzsm)
reclass_displayday <- as.matrix(reclass_displayday)

displaydaypatches<- reclassify(patches,reclass_displayday)

## Plot Spatial Data Output for RZ Storage
#plot(displaydaypatches, xlim = c(280900,282500), ylim = c(4870000,4872000))

{plot(mask(displaydaypatches,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 RZ Soil Moisture Prediction for 1 day", line = -2)
points(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1]), c(ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), pch = 21, col = "black", bg = "red", cex = 2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}

reclass_displayday <- cbind(displayday$patchID,displayday$root.depth)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches2<- reclassify(patches,reclass_displayday)

{plot(mask(displaydaypatches2,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 root.depth Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}

reclass_displayday <- cbind(displayday$patchID,displayday$rz_storage)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches3<- reclassify(patches,reclass_displayday)

{plot(mask(displaydaypatches3,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 rz_storage Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}

reclass_displayday <- cbind(displayday$patchID,displayday$sat_def)
reclass_displayday <- as.matrix(reclass_displayday)
displaydaypatches4<- reclassify(patches,reclass_displayday)

{plot(mask(displaydaypatches4,ws32), xlim = c(275000,276200), ylim = c(3881000,3881600))
title("CW WS32 saturation deficit Prediction for 1 day", line = -2)
scalebar(200, xy=c(275200, 3880800), type='line',divs = "2")}

## Plot Hydrograph
{hydrograph(input=singlerunvalidation$bd, streamflow=singlerunvalidation$bd$observedstreamflow, streamflow2 = singlerunvalidation$bd$streamflow, timeSeries = singlerunvalidation$bd$date, precip = singlerunvalidation$bd$streamflow,
           P.units = "mm", S.units = "mm normalized by basin area", S1.col = 'Blue', S2.col = 'Red')
legend("topleft",inset = .2, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("blue","red"), lty=1:2)}

### end show patch map




patchnumbers<- raster::extract(patches,matrix(c(ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2]), nrow = 3, ncol = 2))


ObserveSoilMoisturePatches <- matrix(c(ObservedSMLocations$Name[3],ObservedSMLocations$Name[4],ObservedSMLocations$Name[5],ObservedSMLocations[[3]][[3]][1],ObservedSMLocations[[3]][[4]][1],ObservedSMLocations[[3]][[5]][1],ObservedSMLocations[[3]][[3]][2],ObservedSMLocations[[3]][[4]][2],ObservedSMLocations[[3]][[5]][2],patchnumbers), nrow = 3, ncol = 4)


## Extract values for a single patch, select patch ID
ObserveSoilMoisturePatches[1,4]
## [1] "136557"
displaypatch<- singlerunvalidation$pd[which(singlerunvalidation$pd$patchID==ObserveSoilMoisturePatches[1,4])]
ObserveSoilMoisturePatches[2,4]
## [1] "141942"
displaypatch2<- singlerunvalidation$pd[which(singlerunvalidation$pd$patchID==ObserveSoilMoisturePatches[2,4])]
ObserveSoilMoisturePatches[3,4]
## [1] "149478"
displaypatch3<- singlerunvalidation$pd[which(singlerunvalidation$pd$patchID==ObserveSoilMoisturePatches[3,4])]

plot(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", main = 'CW WS27 More Detailed RZ Storage Patch #147322', ylim = c(0,.4))

{plot(smws32_1mean$Group.1,smws32_1mean$smois30A, type = "l", ylim = c(0,.4), col = 'darkorange', xlim=as.Date(c("2016-10-1","2019-10-1")), main = "CW WS32 More Detailed Soil Moisture vs Patches",
     xlab = "Date", ylab = "Soil Moisture mm")
lines(smws32_2mean$Group.1,smws32_2mean$smois30A,type = "l", col = 'darkred')
lines(smws32_3mean$Group.1,smws32_3mean$smois30A,type = "l", col = 'blue')
lines(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'red', lty = 2)
lines(displaypatch2$date,(displaypatch2$rz_storage/displaypatch2$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'green', lty = 2)
lines(displaypatch3$date,(displaypatch3$rz_storage/displaypatch3$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", col = 'black', lty = 2)

legend("bottomleft",inset = .01, legend=c("Observed Patch 1","Predicted Patch 1","Observed Patch 2","Predicted Patch 2","Observed Patch 3","Predicted Patch 3"), col=c("darkorange", "red","darkred","green","blue","black"), lty=c(1,2,1,2,1,2))}

{plot(smws32_2mean$Group.1,smws32_2mean$smois30A, type = "l", ylim = c(0,0.4), col = 'BLUE', xlim=as.Date(c("2016-10-1","2019-10-1")), main = "CW WS32 More Detailed Soil Moisture vs Patch #154914",
     xlab = "Date", ylab = "Soil Moisture mm")
lines(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l", xlab = "Date", ylab = "RZ Storage", main = 'CW WS27 More Detailed RZ Storage Patch #154914', col = 'red')
legend("bottomleft",inset = .01, legend=c("Observed Soil Moisture","Predicted Soil Moisture"), col=c("blue","red"), lty=1:1)}

#
#obsws32cal <- obsws32[obsws32$Date >= Caldates[1] & obsws32$Date <= Caldates[2], ]
#obsws32val <- obsws32[obsws32$Date >= Valdates[1] & obsws32$Date <= Valdates[2], ]
#obsws32smcal <- obsws32sm[obsws32sm$Date >= Caldates[1] & obsws32sm$Date <= Caldates[2], ]
#obsws32smval <- obsws32sm[obsws32sm$Date >= Valdates[1] & obsws32sm$Date <= Valdates[2], ]

head(obsws32cal)
##             Date flow discharge_mm observedstreamflow
## 22961 2017-01-01 2.30     2.172645           2.172645
## 22962 2017-01-02 3.04     2.871669           2.871669
## 22963 2017-01-03 2.72     2.569388           2.569388
## 22964 2017-01-04 2.11     1.993165           1.993165
## 22965 2017-01-05 1.92     1.813686           1.813686
## 22966 2017-01-06 1.81     1.709777           1.709777
head(obsws32val)
##             Date flow discharge_mm observedstreamflow
## 23311 2017-12-17 2.14     2.021504           2.021504
## 23312 2017-12-18 2.27     2.144306           2.144306
## 23313 2017-12-19 2.23     2.106521           2.106521
## 23314 2017-12-20 3.39     3.202289           3.202289
## 23315 2017-12-21 2.74     2.588281           2.588281
## 23316 2017-12-22 2.60     2.456033           2.456033
head(obsws32smcal)
##         Date flow discharge_mm observedstreamflow smois30site1a smois30site1b
## 1 2017-01-01 2.30     2.172645           2.172645     0.2194375     0.2468125
## 2 2017-01-02 3.04     2.871669           2.871669     0.2241528     0.2509896
## 3 2017-01-03 2.72     2.569388           2.569388     0.2184132     0.2409236
## 4 2017-01-04 2.11     1.993165           1.993165     0.1952743     0.2236250
## 5 2017-01-05 1.92     1.813686           1.813686     0.1825451     0.2149410
## 6 2017-01-06 1.81     1.709777           1.709777     0.1745729     0.2097500
##   smois60site1a smois60site1b smois30site2a smois30site2b smois60site2a
## 1     0.1988021     0.2463715     0.2162465     0.2141181     0.2191528
## 2     0.2243646     0.2737917     0.2244410     0.2392847     0.2239028
## 3     0.2177813     0.2657361     0.2209722     0.2291771     0.2154028
## 4     0.1993090     0.2483194     0.2058437     0.2039549     0.1967431
## 5     0.1845556     0.2335799     0.1972882     0.1883021     0.1870486
## 6     0.1744132     0.2231181     0.1914861     0.1774236     0.1808646
##   smois60site2b smois30site3a smois30site3b smois60site3a smois60site3b
## 1     0.1775868     0.2830729     0.2573576     0.2625312     0.2399167
## 2     0.2186528     0.2833924     0.2626007     0.2942917     0.2684410
## 3     0.2180486     0.2697604     0.2480833     0.2926944     0.2704132
## 4     0.1968403     0.2399167     0.2206285     0.2794965     0.2609514
## 5     0.1786215     0.2264861     0.2081111     0.2668299     0.2512847
## 6     0.1669340     0.2187014     0.2009896     0.2569028     0.2437292
##   mergedsoilmoisture
## 1          0.2317839
## 2          0.2490255
## 3          0.2422839
## 4          0.2225752
## 5          0.2099661
## 6          0.2015738
head(obsws32smval)
##           Date flow discharge_mm observedstreamflow smois30site1a smois30site1b
## 351 2017-12-17 2.14     2.021504           2.021504     0.1181493     0.2395313
## 352 2017-12-18 2.27     2.144306           2.144306     0.1165625     0.2392708
## 353 2017-12-19 2.23     2.106521           2.106521     0.1135139     0.2332326
## 354 2017-12-20 3.39     3.202289           3.202289     0.1484514     0.2609861
## 355 2017-12-21 2.74     2.588281           2.588281     0.1391528     0.2532083
## 356 2017-12-22 2.60     2.456033           2.456033     0.1312049     0.2428194
##     smois60site1a smois60site1b smois30site2a smois30site2b smois60site2a
## 351     0.1677882     0.2054826     0.1575382     0.1352778            NA
## 352     0.1635625     0.2064132     0.1605903     0.1370000            NA
## 353     0.1600174     0.2057812     0.1585660     0.1363750            NA
## 354     0.1844479     0.2319132     0.1872778     0.1721458            NA
## 355     0.2092639     0.2517917     0.1925139     0.1736111            NA
## 356     0.1954965     0.2400937     0.1847604     0.1649479            NA
##     smois60site2b smois30site3a smois30site3b smois60site3a smois60site3b
## 351    0.09663194     0.2347188     0.1869757     0.2350000     0.2169896
## 352    0.09784375     0.2385694     0.1832847     0.2350000     0.2167847
## 353    0.09828472     0.2289757     0.1792153     0.2350000     0.2160000
## 354    0.12405208     0.2609826     0.2203681     0.2421354     0.2173993
## 355    0.14059028     0.2525278     0.2136979     0.2712326     0.2365937
## 356    0.13011806     0.2397049     0.1995382     0.2683160     0.2404201
##     mergedsoilmoisture
## 351          0.1812803
## 352          0.1813529
## 353          0.1786329
## 354          0.2045600
## 355          0.2121985
## 356          0.2034018
#subset Validation range based on Val Dates set by available SM data

validationsubset <- singlerunvalidation$bd[singlerunvalidation$bd$date >= Valdates[1] & singlerunvalidation$bd$date <= Valdates[2], ]

## Fit tests for streamflow

NSE(validationsubset$streamflow,validationsubset$observedstreamflow)
## [1] 0.151659
NSE(validationsubset$streamflow,validationsubset$observedstreamflow, FUN = log, epsilon = "Pushpalatha2012", na.rm=TRUE)
## [1] 0.6005506
KGE(validationsubset$streamflow,validationsubset$observedstreamflow)
## [1] 0.5055394
# Plot Observed and Predicted Streamflow

{plot(validationsubset$date,validationsubset$streamflow, type = 'l', lty = 2, col = 'RED', main = "Validation Subset Observed and Predicted Streamflow")
lines(validationsubset$date,validationsubset$observedstreamflow, type = 'l', col = 'BLACK')
legend("topright",inset = .1, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:2)}

{plot(singlerunvalidation$bd$date,singlerunvalidation$bd$streamflow, type = 'l', lty = 2, col = 'RED', main = "Observed and Predicted Streamflow")
lines(singlerunvalidation$bd$date,singlerunvalidation$bd$observedstreamflow, type = 'l', col = 'BLACK')
legend("topright",inset = .1, legend=c("Observed Streamflow","Predicted Streamflow"), col=c("black","red"), lty=1:2)}

#Plot Model RZ SM outputs for patch
{plot(displaypatch$date,(displaypatch$rz_storage/displaypatch$root.depth), type = "l",ylim = c(0,0.50), ylab = "Root Zone Storage (mm/m) ", xlab = "Date", main = 'Root Zone Display Patch')
lines(displaypatch$date,(displaypatch$rz_field_capacity/displaypatch$root.depth), type = "l", col = "blue")
lines(displaypatch$date,(displaypatch$rz_wilting_point/displaypatch$root.depth), type = "l", col = "red")
legend("bottomleft",inset = .1, legend=c("Predicted RZ Storage","RZ Field Capacity","RZ Wilting Point"), col=c("black","blue","red"), lty=1:1)}

Model Validation

Valdates
## [1] "2017-12-17" "2018-05-15"
## for now run top run and record values

paramtable[which(KGElist == Reduce(max, KGElist), arr.ind=TRUE),]
##           m        k  soil_dep       m_v      k_v        gw1       gw2
## 65 1.669432 26.16634 0.3261013 0.3147534 11.40159 0.03675462 0.2058475
##          NSE     lnNSE       KGE run     smNSE   smlnNSE      smKGE
## 65 0.4219818 0.3230267 0.6270194  65 -4.489804 -3.079227 -0.1218914
{plot(cwws32_run4$bd$date,cwws32_run4$bd$streamflow, type = "l", col = 'red', main = "Observed Streamflow vs Top KGE Run", xlab = "Date", ylab = "Streamflow")
lines(obsws32$Date,obsws32$discharge_mm, col = 'orange', lwd = 2, lty = 2)}

valsubset <- simmerge4[simmerge4$Date >= Valdates[1] & simmerge4$Date <= Valdates[2], ]

ValKGE<- KGE(sim = valsubset$streamflow, obs = valsubset$observedstreamflow)


valsubsetsm<- simmergesm4[simmergesm4$Date >= Valdates[1] & simmergesm4$Date <= Valdates[2],]

ValsmKGE<- KGE(sim = valsubsetsm$rz_storage/valsubsetsm$rootdepth, obs = valsubsetsm$mergedsoilmoisture)



ValKGE
## [1] 0.5061636
ValsmKGE
## [1] -0.610026